Assessing the impact of plant genetic diversity in shaping the microbial community structure of Vitis vinifera phyllosphere in the Mediterranean

ABSTRACT The aerial surface of the plant (phyllosphere) is the habitat of complex microbial communities and the structure of this microbiome may be dependent on plant genetic factors, local environment or interactions between them. In this study, we explored the microbial diversity present in the phyllosphere of a very diverse set of grapevine cultivars representing the three genetic pools of the species, grown on an experimental plot at Montpellier (French Mediterranean region). We assessed microbiome variation in the phyllosphere using amplicon sequencing of the 16S rRNA gene and of the internal transcribed spacer (ITS), according to the grapevine genetic pools or cultivars, and organs (i.e. leaves and grape berries). The observed microbiome was complex; out of 542 bacterial genera; Pseudomonas, Pantoea, Sphingomonas, and Acinetobacter were the most abundant and almost ubiquitously present across the samples, and out of 267 fungal genera; Aureobasidium, Alternaria, Mycosphaerella and Aspergillus were most represented. Our results illustrated that the microbial taxa were almost uniformly distributed among the genetic pools and only a few cultivar or genetic pool level differences were found, but a very clear differential taxa abundance was found between the leaf and berry samples. Some genus level associations were also observed with certain genetic pools.


Introduction
Vitis vinifera (subsp. vinifera L.), is the main grape species grown for fruit and wine production over the world. It is a natural host of a wide variety of prokaryotic and eukaryotic microorganisms that interact with grapevine, having either beneficial or phytopathogenic effects (Schulz et al. 1999). These microbes also play a major role in fruit yield, grape quality and, ultimately, in the pattern of grape fermentation and wine production (Compant et al. 2011;Bokulich et al. 2016;Belda et al. 2017).
The grapevine phyllosphere is rather less extensively studied as compared to the rhizosphere and endosphere (Vorholt 2012). The phyllosphere (in general) also harbors complex microbial communities involved in many crucial functions such as nitrogen fixation (Jones 1970), carbon sequestration (Bringel and Couée 2015), degradation of pesticides and organic pollutants (Brandl et al. 2001;Kishore et al. CONTACT Prashant Singh singh.prashant.iitkgp@gmail.com All of the data are provided fully in the result section of this paper and the sequence data for 16S and ITS sequences are available at institutional server http://agap-ng6.supagro.inra.fr/inra and can be obtained upon reasonable request to authors. Supplemental data for this article can be accessed here. https://doi.org/10. 1080/21553769.2018.1552628 2005; Bulgarelli et al. 2013). It is a significant and ubiquitous habitat for microorganisms and also an open system that microbes can invade by migration from the atmosphere, soil, other plants and insects (Lugtenberg et al. 2002;Williams et al. 2013). But microbial populations on phyllosphere are also known to live and thrive under harsh environmental factors such as UV radiation, air pollution, temperature fluctuations, water and nutrient availability (Andrews and Harris. 2000;Lugtenberg et al. 2002;Müller and Ruppel 2014). A very fundamental question in microbial ecology is what drivers shape the microbiome on phyllosphere? Environmental conditions at the particular location and biotic factors such as leaf age have been identified as important drivers (Kadivar and Stapleton 2003;Ikeda et al. 2011, Copeland et al. 2015. Some reports on grapevine phyllosphere also suggested that the bacterial and fungal communities of the phyllosphere are minimally affected by the chemical and biological treatments tested, and they mainly differed according to the grapevine location (Gu et al. 2010;Bokulich et al. 2014;Perazzolli et al. 2014). Few authors suggested that in the tropical and temperate forests, the plant genotype is a major driver of the composition of the bacterial communities in the phyllosphere (Lambais et al. 2006;Redford et al. 2010). Another study on Arabidopsis thaliana also illustrated that the plant genetic factors may influence the community composition of the phyllosphere (Bodenhausen et al. 2014).
Until recently, there have been no scientific reports available, analyzing the effect of grapevine genetic factors on the microbiome structure of the phyllosphere. Considering that the microbial diversity present in phyllosphere could be relevant for plant health (Lugtenberg et al. 2002;Compant et al. 2005;Vorholt 2012), a better understanding of how the microbiota associated with grapevine phyllosphere is structured according to the grapevine genetic diversity available at a particular geographic location may provide unexpected opportunities to develop innovative and natural biocontrol methods or phytostimulators against plant pathogen or new breeding scheme for the creation of innovative resistant cultivars. As a first step towards this goal, we explored the bacterial and fungal diversity in the phyllosphere of leaf and berry samples from a set of rather diverse grapevine cultivars that belongs to the three genetic pools of the cultivated grapevine (Nicolas et al. 2016), in the French Mediterranean region. These experiments led us to address two major questions: (i) What microbial diversity is present in the phyllosphere of our Mediterranean vineyards and (ii) how this microbiome structure itself according to the grapevine genetic diversity and plant organs.

Sample collection and DNA extraction
A total of 279 grapevine cultivars were grown in a completely randomized block design at Le Chapitre INRA Villeneneuve-Les-Maguelonne field station near Montpellier (French Mediterranean region). A panel of cultivars representing three genetic pools (western Europe, WW; from eastern Europe, WE; and table grape, TE) was constructed for genome-wide association studies while minimizing relatedness and retaining the main founders of modern cultivated grapevine to optimize the genetic diversity (Nicolas et al. 2016). Nine cultivars were randomly selected from each genetic pool and leaf (with sizes +1 to +4) samples were taken from four to five plants of each cultivar. Leaf samples were taken before spraying of pesticides; each plant had the same age. We collected the leaf samples in the Spring (mid-May 2016) and at the beginning of harvesting season, we also collected samples of berries from the same cultivars. A metadata table containing all the information about the samples and replicates can be downloaded from the GitHub repository (https://github.com/PrashINRA/MetaData_Grapevine Phyllo.git). All samples were washed with an isotonic solution of sodium chloride (0.15 M) containing 0.01% Tween 20 using a horizontal shaker for 1hr at 100 RPM. Afterward, samples were given an ultrasonic bath for 7-10 min using Ultrasonic Cleaner (Branson 5510) for maximum recovery of microbes from the sample surface. The remaining solution was centrifuged at 4,000 g and microbial pellets obtained in a 2-ml Eppendorf tube were collected and stored at −20°C. DNA was extracted from the pellets using the Meta-G-Nome Isolation Kit (Epicentre, Illumina) following the manufacturer's instructions.

PCR amplifications and MiSeQ library preparation.
To access bacterial communities, the V4 region of the 16S ribosomal gene was amplified using primers 515F and 806R (Caporaso et al. 2011). Fungal community diversity and abundance were accessed using modified ITS9 and ITS4 primers targeting the ITS2 region (Blaadid et al. 2013;Lundberg et al. 2013). Two-step PCR was performed to prepare sequencing libraries. PCR1 was designed to perform amplification of the target regions and to add Illumina Nextera transposase sequence to the amplicons. Primers from Illumina kit for dual indexing of the amplicons was used in PCR2. Both forward and reverse primers for PCR1 were amended with frameshift (FS) sequences in their 5 overhang to improve sequence diversity and overall read quality (de Souza et al. 2016). PCR1 was performed in 20 μL reactions with 30 ng of sample DNA using the Advantage 2 PCR kit (Clontech, 639206). PNA PCR clamps were also used to reduce host organelle contamination (de Souza et al. 2016). The same PCR1 was performed for ITS amplification except for the step of PNA annealing. Amplicon replicates were pooled, purified using Agencourt AMPure XP beads (Beckman Coulter) at a bead-to-DNA ratio of 0.7:1, resuspended in 30 μL MilliQ water and evaluated in agarose gels. In PCR2, each cleaned PCR1 product within the same sample received a unique combination of forward and reverse primers (respectively, N7 and S5 Illumina dual index oligos). Afterward, samples were again cleaned using AmPure XP magnetic beads, pooled in equimolar concentrations and sequenced using 2 × 250 bp MiSeq v2 sequencing (Illumina Inc., San Diego, CA, USA).

Data processing and analysis
All RAW data files were imported and processed in the R-environment (R Core Team, 2017) using various codes and inbuilt functions available in different Rpackages. The whole dataset for 16S and ITS amplicon sequences were uploaded and available at the institutional server http://agap-ng6.supagro.inra.fr/inra. Data processing and further analysis were done in two phases. In phase-I, raw data files from both the datasets were filtered and trimmed using the fastq-PairedFilter() function of the dada2 package (Callahan et al. 2016) and bases with low-quality scores were discarded. These filtered files were then processed using the Divisive Amplicon Denoising Algorithm (DADA) pipeline which included the steps of dereplication, core denoising algorithm and merging of the base pairs. Merging function provided global ends-free alignment between paired forward and reverse reads, and merged them together if they overlapped exactly and a table for ribosomal sequence variants (RSVs, a higher analog of operational taxonomic units-OTUs) was constructed, which records the number of times each amplicon sequence variant was observed in each sample. DADA infers sample sequences exactly and resolves differences of as little as one nucleotide (Callahan et al. 2016). Chimeras were removed using the removeBimeraDenovo() function of the dada2 package. OTU sequences were assigned a taxonomy using the RDP classifier and the UNITE database (Wang et al. 2007;Abarenkov et al. 2010) with assignTaxonomy() function of the same dada2 package for 16S and ITS sequences, respectively. Then, at the end of phase-I data processing, a phyloseq data object was created to initiate phase-II data analysis.
In phase-II, a phylogenetic tree for the taxa was constructed using the R-package ape (Paradis et al. 2004) and merged with the phyloseq data object of phase-I. Unassigned taxa and singletons were also removed using the subset_taxa() and prune_taxa() functions of the phyloseq package in R (McMurdie and Holmes 2013). This data object was then used to calculate microbial abundances, α, β diversity analysis and for other statistical tests using various functions in the phyloseq and vegan packages (McMurdie and Holmes 2013; Oksanen et al. 2017).
Prevalence plot for taxa abundances was made using ggplot() function of the ggplot2 package (Wickham 2009) using the entire 16S and ITS data-sets. Chao1 estimates of α diversity (Chao 1987) was measured within sample categories using estimate_richness() function of the phyloseq package. Relative abundances of microbial genera were also plotted using the ggplot2 package (Wickham 2009) on the above data, which were also rarified to even depth of 5,000 reads per sample.
Multidimensional scaling (MDS, also known as principal coordinate analysis; PCoA) was performed using Bray-Curtis dissimilarity matrix (Beals 1984) between samples and visualized by using their base functions in the phyloseq package (McMurdie and Holmes 2013).

Statistical analysis
We analyzed all the data from 16S and ITS amplifications separately in R version 3.3.4 using the dada2, phyloseq and vegan packages. CRAN packages plyr and ggplot2 (Wickham 2009;Wickham 2011) were also used to draw the figures. We assessed the statistical significance (P < 0.05) throughout and whenever necessary, we adjusted P-values for multiple comparisons according to the Benjamini and Hochberg method to control False Discovery Rate (Benjamini and Hochberg 1995) while performing multiple testing on taxa abundance according to sample categories. We performed an analysis of variance (ANOVA) among sample categories while measuring the Chao1 estimates of α-diversity. Stratified permutational multivariate analysis of variance (PERMANOVA) with 999 permutations was conducted on all principal coordinates obtained during PCoA with the adonis() function of the vegan package, to observe the statistical significance of clusters according to the sample categories.
Linear regression (parametric test), and Wilcoxon (Non-parametric) test (Hollander and Wolfe 1973) were performed on taxa abundances against genetic pools using their base functions in R (Myles and Douglas 1973;Bauer 1972).

Quality assessment of the data
Raw demultiplexed sequence data files were generated using high-throughput amplicon sequencing of 16S and ITS ribosomal RNA genes and the number of reads per sample has been taken into account to obtain the depth of the sequencing. Rarefaction curves (number of reads vs number of OTUs) from both the datasets (Figure 1(A and B)) began to level off for most of the samples suggesting a good quality and coverage of both the data-sets and thus we can assume that the microbial communities were reasonably characterized with the sampling effort.

Microbial diversity in the phyllosphere
A total of 5,772,135 16S and 3,807,033 ITS amplicon sequences were generated from 80 samples covering  two sample types (or organ types) and three genetic pools, respectively. We identified 12,875 unique bacterial and 3,413 unique fungal OTUs, in our phyllosphere samples. After removal of unassigned taxa (genus level assignment) and singletons, 6017 unique bacterial and 2075 unique fungal OTUs belonging to 542 bacterial and 267 fungal genera were recovered. Phylum level classification of bacterial and fungal communities was also identified (Figure 2(A and   B)) using the feature prevalence of entire16S and ITS data-sets, which is the number of samples in which a taxon appeared at least once. For example, the phylum Ignavibacteriae had only five unique taxa with the cumulative abundance of thirty-eight and its presence is observed in less than 10% of the samples. Bacterial and fungal communities were heavily dominated by phylum Proteobacteria (relative abundance > 55%) and Ascomycota ( > 65%) respectively.

Effects of genetic diversity on microbial communities in the phyllosphere
Multiple testing on each of bacterial 6017 and fungal 2075 OTUs (with adjusted p-values:adjp and controlling false discovery rates) was performed according to cultivars and genetic pools and apart from two bacterial taxon (OTU1309, genus: Gemmatimonas, adjp = 0.0209, FDR = 0.06017 and OTU120, genus: Hymenobacter, adjp = 0.036, FDR = 0.05), and one fungal taxon (OTU63, genus: Penicillium, adjp = 0.02, FDR = 0.028), which were differentially abundant between WW, and WE, we did not recover any taxa whose abundance is significantly different (statistically) among the genetic pools.
Relative abundances for the twenty most abundant genera were plotted as well for each cultivar within their genetic pools (Figure 3(A and B)) and microbial genera were quite uniformly abundant among the three genetic pools. This pattern was also the same when we analyzed the abundances in leaf and berry samples individually within these three genetic pools (Figure 4(A and B)), except for few cultivar level differences (e.g. bacterial genus Vagococcus in the cultivars of TE and fungal genus Pichia in the cultivars of WW genetic pool). To test the association of these genera with genetic pools, we performed a linear regression for abundances of these genera against genetic pools (parametric test). As the Pichia also seems more abundant in the phyllosphere of berries (Figure 4(B)), we also added this as confounders to the regression and the results indicated a highly significant association of these genera to TE and WW genetic pools, respectively (Tables 1 and 2). We also observed that the abundance data for these genera were not normally distributed and therefore performed a nonparametric test (Wilcoxon rank sum test), that confirmed the association.
The Chao1 estimator of alpha diversity was also measured and plotted according to the genetic pools (Figure 3(C and D)) and again we did not observe a very significant genetic pool wise differences in these estimates (ANOVA, for 16S data: Chao1, P = 0.033; for ITS data: Chao1, P = 0.041).
Microbial community structure assemblages among the three genetic pools were also compared using PCoA to look for the genetic pool wise patterns of microbiota present in the phyllosphere. Taxa in both the PCoA plot (Figure 3(E and F)) were clustered together (PERMANOVA, for 16S data: at F = 0.971, R 2 = 0.285, P = 0.408; for ITS data: at F = 0.991, R 2 = 0.172, P = 0.394), which also indicated the impact of genetic diversity is less evident. Results were the same when PCoA was performed on the data-sets grouped within 27 grapevine cultivars (Supplementary data S1).

Effect of organs on phyllosphere microbiome
Multiple testing on taxa abundances in the phyllosphere of leaves and berries gave 17 bacterial and 33 fungal OTUs whose abundance was significantly different between these two organs. The data revealed the organ-specific patterns of phyllosphere microbiota in these grapevine cultivars. Tables 3 and 4 are provided for 16S and ITS data respectively to display various bacterial and fungal OTUs (along with their respective genera) with their false discovery rates (FDRs) and adjusted p-values. According to the corrected p-values and FDRs, 5 bacterial (e.g. Pseudomonas and Pantoea. adjusted P-value; adjp = 0.0038 & FDR = 0.00118) and 31 fungal genera (e.g. Aspergillus and Mycosphaerella. adjp = 0.0005 & FDR = 0.000129) were most differed between leaf and berries were, respectively.
Relative microbial abundances for top twenty taxa was also calculated on leaf and berry samples (also grouped in genetic pools; Figure 4(A and B)) and differential abundances on both sample type were clearly visible. Leaf phyllosphere was heavily occupied by bacterial and fungal genera of Pseudomonas and Pantoea & Aureobasidium, Mycosphaerella, respectively. On the other hand, berry surfaces mainly comprised of bacterial genera of Acinetobacter and Sphingomonas & with fungal genera of Aureobasidium, Aspergillus and Pichia.
To investigate the influence of leaves and berries, we also compared Chao1 estimates of alpha diversity between leaf and berry samples (Figure 4(C and D)) and these estimates were also significantly different (ANOVA, for 16S data: Chao1, P = 0.007; for ITS data: Chao1, P = 4.53e-08).

Figure 3.
Relative abundances of (A) bacterial and (B) fungal genera present on each cultivar, grouped within their genetic pools (9 cultivars per genetic pool, top 20 taxa, characterized to the genus level and datasets were rarified to 5000 sequence reads per sample). Chao1 estimates of α-diversity for (C) bacterial and (D) fungal data-sets for three genetic pools. PCoA plots using Bray-Curtis distance between samples for (E) bacterial and (F) fungal data-sets among three genetic pools, explaining > 60% variations with first two axes (taxa with variance < 1e-05 were trimmed).

Discussion
Our analysis based on high throughput 16S and ITS profiling identifies the presence of complex microbiota in the phyllosphere of leaves and fruits (berries) of grapevine cultivars grown in our Mediterranean vineyard and it is dominated by bacterial genera of Pseudomonas, Sphingomonas, Enterobacter and the fungal genera of Aureobasidium, Alternaria, Cladosporium, respectively which is concordant with the findings of other grapevine related studies (Zarraonaindia et al. 2015;Zhang et al. 2017). High relative abundances of some other microbial genera such as Pantoea and Mycosphaerella have also been identified in this study, which has been reported in few grapevine cultivars before as endophytes (Bell et al. 1995;Baldan et al. 2015). This is not uncommon as epiphytes and endophytes are separated by a thin boundary between their habitats and due to vertical and horizontal microbial . Relative abundances of (A) bacterial and (B) fungal genera present on leaf and berry samples, also grouped within their genetic pools (top 20 taxa, characterized to the genus level, datasets were rarified to 5000 sequence reads per sample). Chao1 estimates of αdiversity for (C) bacterial and (D) fungal data-sets for both the organ types. PCoA plots using Bray-Curtis distance between samples for (E) bacterial and (F) fungal data-sets as per leaf and berry samples based on Bray-Curtis distance matrices, explaining > 60% variations with first two axes (taxa with variance < 1e-05 were trimmed).
transfers (Frank et al. 2017) sharing of the major chunk of OTUs are inevitable (Bodenhausen et al. 2013). The bacterial genus Vagococcus has also not been widely reported in plants by research communities except rhizosphere and phyllosphere of Rice (Mwajita et al. 2013). Hence, the specific abundance of this genus must be further identified to have the preliminary view of its functionality in the Mediterranean vineyards. The reason for its abundance could be the interaction between plant genetic factors and environmental conditions at this specific geographic location of the vineyard. Epiphytes (the phyllosphere microbes) associated with grapevine have been suggested to originate from soil, but are distinct from those in the rhizosphere microbiome (Zarraonaindia et al. 2015); this is most likely a consequence of the physio-chemical composition and surrounding environment which strongly modulates microbial  community structure and its dynamics (Wagner et al. 2016). Other reports also evidenced that both environment and the plant genotype could be the major drivers for epiphytic community structuring (Redford et al. 2010;Turner et al. 2013). Our preliminary study based on random sampling of 27 cultivars (from three genetic pools) also indicated that there is probably an impact of grapevine genetic diversity over the microbial composition in the phyllosphere, but it is not quite evident as we found only a few microbial OTUs were differentially abundant among genetic pools. Sampling from cultivars (among these genetic pools) which are more distant in the context of their genetic relatedness must be done in the future to further explore the impact of this genetic diversity. Few genera were also found specifically associated with some cultivars of certain genetic pools (e.g. Vagococcus with the genetic pool TE) and this association (if further confirmed with above mentioned sampling strategies), should be taken into account in developing new selective breeding strategies in order to have the putative beneficial role of the phyllosphere as a performance trait of the cultivars. Moreover, environmental control in shaping microbiome in these genetically diverse grapevine cultivars should also be further investigated by sampling at different geographic locations displaying variable climatic conditions.
On the other hand, leaf and berry samples clearly displayed very distinct microbial patterns. Type of genera present and their taxa abundances was significantly different in both the organs. The physical features of berry surface like the number of waxy layers (or bloom, which prevent water loss from the skin) and their thicknesses are cultivar-specific (Knoche and Lang 2017). These physical features could influence the contact and permeability of the grape berry cuticle to different microorganisms as observed for some pathogens, such as Botrytis cinerea (Herzog et al. 2015) and could be the reason for organ-specific microbiome differences and deserve further investigations.
This finding is also consistent with the few other findings, in which organ-specific microbial patterns have been reported in sugarcane (de Souza et al. 2016) and in some commercially important grapevine cultivars (Bokulich et al. 2014). Our leaf samples were majorly occupied by the bacterial and fungal genus Pseudomonas, Pantoea and Sphingomonas & Aureobasidium, Mycosphaerella and Cladosporium. At the other end, berry surfaces displayed higher abundances of bacterial genus like Acinetobacter, Gluconobacter, Enterobacter, but major fungal abundances were similar in both leaf and fruit surfaces except the genus of Aspergillus and Pichia. Pichia (a yeast, family Saccharomycetaceae) was also found specifically abundant in berries of grapevine cultivars of the genetic pool WW. Taxonomy of Pichia is not fully resolved, and thus, a large diversity of roles in winemaking may be expected within this genus with some species inducing potential faults in winemaking (Fugelsang and Edwards 2010). Therefore the information regarding its association with certain genotypes should further be investigated in the context of wine fermentation.
A richly diverse fungal component of the grapevine microbiome has also been uncovered in this work and it could also be particularly significant because there is not sufficient information on the potential risks or benefits of plant-fungi associations. Grapevine associated microbial communities are relevant to industrial fermentation processes for wine production. Based on our results it can be assumed that grape juice used for wine production harbors a diverse bacterial and fungal community originating from its phyllosphere as well (e.g. Pantoea and Aspergillus).
Some species of most abundant microbial genera we found (e.g. Pseudomonas and Mycosphaerella) have been previously reported for acting as biocontrol agents or BCAs (Kurose et al. 2016;Jousset et al. 2006). An interesting question would be to evaluate how to integrate microbial community studies into traditional biocontrol approaches? This integration could provide a better understanding of how microbial communities are interacting with each other, with the host plant, pathogen or with BCAs, which would be definitely helpful for designing a novel biocontrol method. This also suggests that there is an open field for further studies of the possible role of bacterial and fungal colonizers in plant growth, development and response to biotic and abiotic stress.
A whole-genome shotgun sequencing followed by metagenomic analysis (Qin et al 2010) can add a more detailed layer of information to the taxonomical characterization of a wide variety of grapevine samples, by generating information on the gene composition of the bacteria and fungi present. This information can, in turn, be used to discover new genes and to formulate putative functional pathways and modules, thus could provide insight into functional and genetic microbiome variability. Apart from metagenomics, the use of additional tools such as RNA-Seq (for meta-transcriptomics) may offer a more informative perspective as it can reveal details about populations that are transcriptionally active and not just identify the taxa and genetic content of microbial populations. Moreover, the integration of different omic approaches (e.g. meta-transcriptomics & meta-proteomics) may open a window into discovering the regulatory mechanisms orchestrating observed gene expressions, thereby uncovering how host-microbe and microbe-microbe interactions that regulate microbiome activity. research; Prashant Singh and Sylvain Santoni performed the lab experiments; Prashant Singh and Alex Gobbi analyzed the data; Prashant Singh, Jean-Pierre Péros and Patrice wrote this paper.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This study was funded by the Horizon 2020 programme of the European Commission within the Marie Skłodowska-Curie Innovative Training Network 'MicroWine' [grant number 643063].