Metagenomic Characterization Reveals Pronounced Seasonality in the Diversity and Structure of the Phyllosphere Bacterial Community in a Mediterranean Ecosystem

We explore how the phyllosphere microbial community responds to a very seasonal environment such as the Mediterranean. For this, we studied the epiphytic bacterial community of a Mediterranean ecosystem in summer and winter, expecting to detect seasonal differences at their maximum. With high-throughput sequencing (HTS), we detected the operational taxonomic units (OTUs) present in the phyllosphere and also in the surrounding air. The epiphytic community is approximately five orders of magnitude denser than the airborne one and is made almost exclusively by habitat specialists. The two communities differ considerably but Proteobacteria and Actinobacteria are dominant in both. Of the five most abundant phyllosphere OTUs, two were closely related to Sphingomonas strains, one to Methylobacterium and the other two to Rhizobiales and Burkholderiales. We found the epiphytic community to become much richer, more distinct, even and diverse, denser and more connected in summer. In contrast, there was no difference in the level of bacterial colonization of the phyllosphere between the two seasons, although there were seasonal differences for individual taxonomic groups: Firmicutes, Gemmatimonadetes and Chlroroflexi had a higher participation in summer, whereas the major Proteobacteria classes presented reverse patterns, with Betaproteobacteria increasing in summer at the expense of the prominent Alphaproteobacteria.


Introduction
Life on leaf surfaces that are exposed to extreme variations in meteorological and other environmental factors looks like the life of a fugitive. However, existing evidence suggests that microbial communities are well established, that only a subset of the air-introduced taxa can colonize them and that leaves are not passive acceptors of microbes deposited from the air [1][2][3]. There is also air community. These explorations allow a better understanding of the structure of the microbial community of the Mediterranean phyllosphere, of the dominant players and their strategies, and also of their responses to the marked seasonality of the Mediterranean environment.

Study Site
The study site is in close proximity to the sea, at Halkidiki, in northern Greece (40 • 09 N, 23 • 54 E), an area with a Mediterranean-type climate; July is the hottest month of the year. To assess the phyllosphere and airborne microbial communities, we sampled from an area of less than 5 ha, in summer (July) and winter (January). In both cases, sampling started approximately 1.5 h after dawn on rainless and practically windless days and lasted for approximately 3 hrs. We collected weather data for the days of sampling from the weather station of Neos Marmaras, which is the nearest station to our sampling site, at a distance of 13 km on a straight line. We also took measurements of temperature and wind speed on the spot, at the end of sampling (Supplement Table S1).

Sampling and Sample Processing
To assess the epiphytic microbial community of this Mediterranean ecosystem, we took leaf samples from nine perennial plant species that differ in plant habit and associated leaf traits [35] and are representative of the local vegetation and the phyllosphere habitats that are available for microbial growth. These are the evergreen sclerophyllous Arbutus unedo, Myrtus communis, Phillyrea latifolia, Pistacia lentiscus and Quercus coccifera, the seasonal dimorphic Cistus incanus and Lavandula stoechas, and the herbaceous Calamintha nepeta and Melissa officinalis. The sampled area consisted of a flat land, at the mouth of a small torrent that is dry in summer, and of two low-hill slopes at each side: the more mesic and densely covered NE facing slope, dominated by evergreen sclerophyllous species, and the more xeric SW slope, dominated by seasonally dimorphic species. The herbaceous species were collected from the flat land.
For each selected species, at random, we collected mature leaves from five mature individuals. Approximately 0.3 g of the leaves of each sampled individual was weighed and put in 1.2 mL of sterile phosphate buffer (1XPBS: 137 mM NaCl, 10 mM phosphate, 2.7 mM KCl, pH 7.4). Samples were then sonicated in an ultrasonic cleaner for 10 min with the temperature of water not exceeding 20 • C. Buffer suspensions were centrifuged (3-30K Sigma GmbH, Osterode am Harz, Germany) at 9500 rpm, at 4 • C, for 20 min. Pellets of the same species were kept at −20 • C and analyzed as one sample in DNA extraction.
Air sampling took place at the flat land, close to the populations of the herbaceous species, at almost equal distance from the populations of the evergreen schelophyllous and the seasonally dimorphic species that were sampled. For this, we used an Andersen six-stage microbial impaction sampler (Andersen 2000 Inc., Atlanta, GA, USA) deployed on a tripod at a height of 1.5 m above ground. Air sampling was carried out twice for each sampling day, once in the beginning and once at the end of the phyllosphere sampling. Airborne bacteria were collected by means of vacuum filtration for 15 min at a flow rate of 28.3 L min −1 on sterile cellulose filters (0.22 µm pore size, 90 mm diameter, MF-Millipore; Merck KGaA, Darmstadt, Germany) that were placed on top of the Andersen plates. These filters were then aseptically put into test tubes containing 10 mL sterile phosphate buffer (same as for the leaf samples) and transferred in an icebox to the Lab. They were sonicated in an ultrasonic cleaner (Transsonic 460 Elma, Partelli, Brescia, Italy) for 10 min, as described above, vortexed for 2 min and then the filters were removed. Buffer suspensions were centrifuged as above. Pellets were kept at −20 • C and were analysed as one sample per sampling day in DNA extraction.

DNA Extraction, Composition and Abundance of Microbial Communities
Total genomic DNA was extracted by using the PowerSoil DNA isolation kit from Qiagen Laboratories (Hilden, Germany) in accordance with the manufacturer's instructions. The quantity of the DNA was between 5 and 10 ng µL −1 , as measured by Nanodrop (ND-1000, Thermo Fisher Scientific, Waltham, MA, USA). The DNA samples were amplified using the bacterial-specific primers 785F (5 -GGATTAGATACCGTGGTA-3 ) and 1185mR (5 -GAYTTGACGTCATCCM-3 ), which target and amplify a final domain of around 415 bp of the 16S rRNA gene in theV4-V6 SSU region [36,37]. The PCR conditions were as follows: 94 • C for 3 min followed by 28 cycles at 94 • C for 30 s, 53 • C for 40 s and 72 • C for 1 min, which was followed by a final elongation step at 72 • C for 5 min. After amplification, PCR products were checked in 2% agarose gel to determine the success of amplification and the relative intensity of the bands. Subsequently, the PCR products were purified using calibrated Ampure XP beads (Beckman Coulter Life sciences, Brea, CA, USA) and the purified products were used to prepare the DNA libraries by following the Illumina MiSeq DNA high-throughput library preparation protocol. DNA library preparation and sequencing were performed at Mr. DNA (www.mrdnalab.com; Shallowater, TX, USA) on a MiSeq following the manufacturer's guidelines.
For the quantification of bacterial 16S rRNA genes in our samples, a primer set was used as previously described [38]. A serial of 10-fold dilutions of a recombinant plasmid containing a partial fragment of a bacterial 16S rRNA gene was used as external standard to perform the standard curve. The standard dilutions ranged from 10 4 to 10 10 . The real-time PCR was performed in a LightCycler 480 (Roche Basel, Switzerland) instrument using the LightCycler 480 SYBR Green Master I (Roche, Basel, Switzerland)) following the manufacturer's instructions. All samples, standards and negative controls were tested in triplicates. Finally, we used cycle threshold (CT) values to determine the number of 16S rRNA gene copies in our samples.

Read Processing
The produced reads were processed using the mothur v1.34.0 software [39], following the standard operating procedure [40]. Briefly, forward and reverse reads were joined, and the barcodes were removed. Reads below 200 bp, with homopolymers higher than 8 bp and with ambiguous base calls, were removed from downstream analysis. The remaining reads were de-replicated to the unique sequences and aligned independently against SILVA 128 database, containing 1,719,541 bacterial SSU rRNA sequences [41]. Then, the reads suspected for being chimeras were removed using the UCHIME software [42]. Almost 25% of the reads were removed at this step. The remaining reads were clustered into operational taxonomic units (OTUs) at 97% similarity. To obtain a rigorous dataset, OTUs with a single read in the entire dataset were removed from the analysis as they were suspected of being erroneous sequences [43]; these corresponded to almost 50% of the OTUs. The resulting dataset was normalized to the lowest number of reads (5038 reads) with the subsample command in mothur. Taxonomic classification was assigned using SINA searches on the Silva 128 curated database [44]. After normalization, the reads belonging to OTUs affiliated with unclassified sequences at the domain level were removed in order to be confident that the produced dataset includes only bacterial reads. Chloroplast-related reads that were recovered were also removed from the dataset; there were no mitochondrial OTUs. Raw reads were submitted to GenBank-SRA under the accession number SRX3316528.

Data Analysis
We used the normalized dataset for all our analyses. Shannon, Simpson and Pielou's evenness alpha-diversity estimators were calculated with the PAST 2.17c software [45]. The bacterial assemblages of the different samplings were compared using the Plymouth routines in the multivariate ecological research software package PRIMER v.6 [46]. The Bray-Curtis dissimilarity coefficients were calculated to develop the matrix based on OTU abundance in order to identify interrelationships between the samples and construct the cluster plots. For the determination of the OTUs responsible for the within and between group dissimilarities, the similarity percentage analysis (SIMPER) was applied [47].
We used paired t-tests to compare the summer and winter values of (i) the diversity indices of the epiphytic bacterial community, (ii) the richness (number of different OTUs belonging to bacteria) and abundance (bacterial 16S rRNA gene copies g −1 plant tissue) of the entire epiphytic bacterial community, and (iii) the richness and abundance of OTUs corresponding to the bacterial phyla and to the classes of the dominant Proteobacteria that were represented in the epiphytic community. For this analysis, we used Statistica 7 for Windows (StatSoft, Tulsa, USA).
OTUs were classified as abundant or rare in relation to their overall relative abundance, following HTS studies on prokaryotes [48][49][50]. Abundant OTUs were defined as those with relative abundance > 1%, whereas rare OTUs as those with abundance < 0.1% of the total number of reads in the entire dataset. As this classification operates at two levels, OTUs could be also distinguished as locally rare or locally abundant after the total number of reads per sample.
Furthermore, OTUs were classified as generalists and specialists based on Levins' index [51]. Levins' proposed that niche breadth could be estimated by measuring the individuals' uniformity of distribution among the resource states. For this, specialization of each individual OTU was calculated according to Pandit et al. [52], using Levins' niche width (B j ) index [51] (1): where pij is the proportion of OTU j in sample i, and N is the total number of samples. Therefore, B j describes the extent of niche specialization based on the distribution of OTU abundances without taking into account the abiotic conditions in a local community. The values of the index range between 1 for singletons and a maximum value that varies depending on the dataset, which in our case was 15 (top generalist). OTUs with B j index higher than 10 were arbitrarily considered as generalists, while OTUs with B j lower than 5 as specialists [53].
For the network analysis, the relationship between OTUs was characterized through MINE statistics by computing the Maximal Information Coefficient (MIC) between each pair of OTUs, based on the number of reads of each OTU [54]. MIC captures correlations between data and provides a score between 0 (no correlation) and 1 (very strong correlation) that represents the strength of a relationship between data pairs. The sample pairings with MIC values > 0.5 (corresponding to p-values < 0.05) were used to visualize the strong correlations between OTUs [54]. The topological parameters of the respective networks were calculated with Cytoscape 3.0 [55]. Network Randomizer 1.1.2 [56] was used to generate random networks of the dataset in order to compare with and assess air and phyllosphere networks' density.

Results
As determined by quantitative real time PCR, the number of bacterial 16S rRNA gene copies per gram of plant tissue of the Mediterranean phyllosphere studied was of the order of magnitude of 10 8 and did not differ between seasons ( Table 1). The number of bacterial 16S rRNA gene copies per cubic meter of air was of the order of 10 6 . Values for microbial abundance in the air and in the phyllosphere cannot be directly compared as the units differ-i.e., in number of copies per unit of weight, in the case of the phyllosphere, and per unit of volume, in the case of the air. Nevertheless, taking into consideration that the density of air at sea level is approximately 1/800th the density of water, what corresponds to 1.25 kg m −3 [57], we can estimate the microbial abundance in the air on a per unit of weight basis (Table 1). Values are of the order of 10 3 per gram of air.
After read processing, the removal of singletons, non-bacterial sequences and normalization, a total of 890 different OTUs were detected. Rarefaction curves that were calculated for all samples (Supplement Figure S1) indicate that a satisfying part of the diversity was recovered with the sequencing effort applied in most samples. There were 771 OTUs in the overall summer dataset and 430 in the winter one. Of these, 750 and 420 OTUs, respectively, corresponded to the phyllosphere. Air samples had far lower numbers of OTUs: 128 in summer and 86 in winter. The Shannon diversity index of the epiphytic community, estimated after the OTU abundances, was higher in summer (3.63 ± 0.18) than in winter (2.95 ± 0.15). Similarly, the Pielou's evenness index was higher in summer (0.68 ± 0.02 vs. 0.59 ± 0.03), but the Simpson diversity index did not differ between seasons (Table 1). On average, there were 233.8 ± 22.4 OTUs per phyllosphere sample in summer, whereas 146.0 ± 10.9 in winter (Table 1). Summary results regarding the different OTUs in leaf and air samples are given in Table 2.
Proteobacteria was the prominent phylum in the bacterial community of the Mediterranean phyllosphere examined (Table 3), followed by Actinobacteria. Bacteroidetes, Firmicutes, Acidobacteria, Deinococcus, Chloroflexi and Planctomycetes had a considerable participation, though two to three orders of magnitude lower than that of the previous two phyla. Other participating phyla were Gemmatimonadetes, Armatimonadetes, Verrucomicrobia and also Latescibacteria, Nitrospirae, Saccharibacteria, Spirochaetae and the WD272 group; there were also a few unidentified bacteria. Within Proteobacteria, Alphaproteobacteria was the dominant class, followed by Betaproteobacteria. The rank at the phylum and class levels remained quite similar if instead of OTU abundance, the number of different OTUs representing each of the above taxa (OTU richness) was taken into consideration ( Figure 1).    At a lower taxonomic level, Rhizobiales was the bacterial order represented in the phyllosphere by the highest number of OTUs (89) ( Table 4). Rhodospirillales, Micrococcales, Sphingomonadales, and Burkholderiales followed, all with similar number of different OTUs (47)(48)(49)(50)(51)(52)(53)(54)(55)(56)(57). Rhizobiales was first in rank in the air, too, but represented by a far lower number of OTUs (18); again Burkholderiales, Sphingomonadales and Micrococcales followed, represented by [12][13][14][15][16] OTUs. There were also striking differences. Representatives of Acidimicrobiales and Planktomycetales, with a high participation in the phyllosphere (18 and 16 OTUs, respectively), were absent in the air, whereas other major taxa were at a far lower rank in the air than in the phyllosphere and vice versa ( Table 4).
The most abundant OTU in the phyllosphere was closely related to a Sphingomonas strain. Of the next four most abundant OTUs, one was closely related to another Sphingomonas strain, one to a Methylobacterium strain and the other two to Rhizobiales and Burkholderiales strains. At a lower taxonomic level, Rhizobiales was the bacterial order represented in the phyllosphere by the highest number of OTUs (89) ( Table 4). Rhodospirillales, Micrococcales, Sphingomonadales, and Burkholderiales followed, all with similar number of different OTUs (47)(48)(49)(50)(51)(52)(53)(54)(55)(56)(57). Rhizobiales was first in rank in the air, too, but represented by a far lower number of OTUs (18); again Burkholderiales, Sphingomonadales and Micrococcales followed, represented by [12][13][14][15][16] OTUs. There were also striking differences. Representatives of Acidimicrobiales and Planktomycetales, with a high participation in the phyllosphere (18 and 16 OTUs, respectively), were absent in the air, whereas other major taxa were at a far lower rank in the air than in the phyllosphere and vice versa (Table 4). Table 4. The first fifteen orders of phyllosphere bacteria ranked in descending order by the number of OTUs detected in the phyllosphere of the Mediterranean ecosystem studied that belong to them. Given is also the representation of these orders in the surrounding air. Empty spaces under 'Air' mean that taxa other than the ones presented here had higher representation in the air. The most abundant OTU in the phyllosphere was closely related to a Sphingomonas strain. Of the next four most abundant OTUs, one was closely related to another Sphingomonas strain, one to a Methylobacterium strain and the other two to Rhizobiales and Burkholderiales strains.

Orders
The vast majority of OTUs (89%) in the whole database were categorized as rare. Only 15 OTUs were categorized as abundant (Tables 2 and 5); these were present in all phyllosphere habitats and in the air and were all locally abundant in one and up to four habitats. There were 82 universal OTUs that were detected in all phyllosphere habitats ( Figure 2). Corresponding to 9% of all recovered OTUs, these universal OTUs belong to Proteobacteria, Actinobacteria, Bacteroidetes, Firmicutes and Acidobacteria. Markedly higher was the number of OTUs (381) that were found in only one phyllosphere habitat ( Figure 2). These narrow-niche OTUs made 44% of the phyllosphere total (Table 2).
According to the Levins' index, the vast majority of OTUs (89%) were specialists, with only 1% being generalists. All 10 generalist OTUs belong to Alphaproteobacteria and Actinobacteria; these were also among the 82 universal OTUs and, with the exception of one, they were also present in the air (Table 5). Of the specialist OTUs, six were abundant. None of the members of Alphaproteobacteria, a major bacterial class in the Mediterranean epiphytic community both in terms of OTU richness and OTU abundance, participated in the group of abundant specialists.   According to the Levins' index, the vast majority of OTUs (89%) were specialists, with only 1% 8 being generalists. All 10 generalist OTUs belong to Alphaproteobacteria and Actinobacteria ; these 9 were also among the 82 universal OTUs and, with the exception of one, they were also present in the air (Table 5). Of the specialist OTUs, six were abundant. None of the members of 11 Alphaproteobacteria, a major bacterial class in the Mediterranean epiphytic community both in terms 12 of OTU richness and OTU abundance, participated in the group of abundant specialists.

13
There was a pronounced seasonal difference in the structure of the phyllosphere microbial 14 community. OTU richness was higher in summer for the entire community (Table 1), for all major 15 bacterial phyla, i.e., Proteobacteria Actinobacteria, Bacteroidetes, Firmicutes, Acidobacteria, and seasons ( Figure 1). For the taxa differing between seasons in terms of OTU richness, the summer

19
value was in general one to three times higher than the winter value, but Chloroflexi was essentially 20 a summer taxon; all but one of the Chloroflexi OTUs were detected in summer. Gemmatimonadetes 21 was similarly a summer taxon but it was represented by far fewer OTUs than Chloroflexi. In terms 22 of OTU abundance, there was no difference between seasons of the phyla represented, except for

23
Chloroflexi, Gemmatimonadetes and Firmicutes that were more abundant in summer. Also, the class 24 of Betaproteobacteria was more abundant in summer, whereas that of Alphaproteobacteria in winter 25 (Table 3). In summer, there was also a higher number of OTUs with exclusive occurrence in only one 26 phyllosphere habitat (Table 1).

27
Given that our sampling did not allow direct comparisons between the phyllosphere and the There was a pronounced seasonal difference in the structure of the phyllosphere microbial community. OTU richness was higher in summer for the entire community (Table 1), for all major bacterial phyla, i.e., Proteobacteria Actinobacteria, Bacteroidetes, Firmicutes, Acidobacteria, and Chloroflexi, for other minor ones, such as Deinococcus and Gemmatimonadetes, and also for the Proteobacteria classes, except for the dominant Alphaproteobacteria that did not differ between seasons (Figure 1). For the taxa differing between seasons in terms of OTU richness, the summer value was in general one to three times higher than the winter value, but Chloroflexi was essentially a summer taxon; all but one of the Chloroflexi OTUs were detected in summer. Gemmatimonadetes was similarly a summer taxon but it was represented by far fewer OTUs than Chloroflexi. In terms of OTU abundance, there was no difference between seasons of the phyla represented, except for Chloroflexi, Gemmatimonadetes and Firmicutes that were more abundant in summer. Also, the class of Betaproteobacteria was more abundant in summer, whereas that of Alphaproteobacteria in winter (Table 3). In summer, there was also a higher number of OTUs with exclusive occurrence in only one phyllosphere habitat ( Table 1).
Given that our sampling did not allow direct comparisons between the phyllosphere and the airborne microbial communities, we subjected our samples to Cluster analysis, according to Bray-Curtis dissimilarities. Three major clusters emerged at >20% level of similarity (Figure 3). Cluster A consisted of the two air samples. Clusters B and C, separated at a 35% level of similarity, consisted of phyllosphere samples, primarily of winter and summer, respectively. The relative participation of the phyla and Proteobacteria classes in each of these three clusters is shown in Figure 4. Bacteria that do not belong to Proteobacteria, Actinobacteria, Bacteroidetes and Firmicutes had a very low contribution in the air, less than 1%, and whereas Alphaproteobacteria was the dominant Proteobacteria class in the phyllosphere, Betaproteobacteria was the dominant one in the air. Also, Actinobacteria had a far higher relative participation in the summer phyllosphere samples.

43
The SIMPER analysis, based on within cluster similarities, indicated the assemblages of OTUs 44 that were responsible for the formation of each cluster: 17 OTUs for cluster A, 24 for cluster B, and 44 45 for cluster C (Supplement Table S2). Seven OTUs were responsible solely for the formation of cluster 46 A. Of these, three belonged to Alphaproteobacteria, two to Firmicutes and the remaining two to unknown strains that had been reported to be isolated from the soil environment. For the formation 48 of cluster B, three OTUs were solely responsible, of which two belonged to Alphaproteobacteria and 49 one to Acidobacteria. For the formation of cluster C, 19 OTUs were solely responsible, with > 50% of The SIMPER analysis, based on within cluster similarities, indicated the assemblages of OTUs that were responsible for the formation of each cluster: 17 OTUs for cluster A, 24 for cluster B, and 44 for cluster C (Supplement Table S2). Seven OTUs were responsible solely for the formation of cluster A. Of these, three belonged to Alphaproteobacteria, two to Firmicutes and the remaining two to unknown strains that had been reported to be isolated from the soil environment. For the formation of cluster B, three OTUs were solely responsible, of which two belonged to Alphaproteobacteria and one to Acidobacteria. For the formation of cluster C, 19 OTUs were solely responsible, with > 50% of them belonging to Actinobacteria strains (Supplement Table S2).
The correlations of the OTUs responsible for the formation of each cluster were calculated according to Maximal Information Coefficient (MIC) values. Network analysis was implemented on the OTUs exhibiting strong positive correlations (high MIC values, corresponding to a p-value < 0.05) in each cluster ( Figure 5). Even though there are no physical interactions among OTUs found on the leaves of different host plants, network analysis provides indications of similar or opposing trends in OTU occurrence patterns. From the topological parameters calculated for each network, it was evident that the microbial community associated with cluster C was the most connected of the three.
In particular, the clustering coefficient, the number of shortest paths, and the average number of neighbors took the highest values in cluster C (Table 6).

64
In particular, the clustering coefficient, the number of shortest paths, and the average number of neighbors took the highest values in cluster C (Table 6). Network analysis also provided a visualization of the OTU correlations, indicating that different taxonomic groups had strong correlations both within and among them in each cluster. In cluster B, consisting mainly of winter samples, Alphaproteobacteria showed the highest number of correlations followed by Actinobacteria ( Figure 5). In cluster C, consisting mainly of summer samples, Actinobacteria presented the highest number of correlations, while Betaproteobacteria, with a negligible role in the network structure for the other two clusters, presented high connectivity here. The topological parameters of the networks suggest that cluster C has a denser network structure than the other clusters: the Clustering Coefficient and the Average Number of Neighbors were higher in the network of cluster C (Table 6). Furthermore, topological parameters, such as the Clustering Coefficient, Centralization, and Heterogeneity, were approximately 2.5, 2 and 1.8 times higher, respectively, for the observed cluster C network relative to a random network of the same size suggesting indeed a denser structure than expected by chance alone [58].   Table S2. The size of the nodes is analogous to the clustering coefficient of each OTU with larger nodes representing key OTUs for the network, in terms of connectivity and centrality.

Discussion
The study of the phyllosphere microbial community in a Mediterranean ecosystem showed that Proteobacteria is the dominant phylum in terms of both richness and abundance and that this holds true irrespective of season. This is also true for the airborne microbial community. However, the relative abundance in terms of average number of reads of bacterial classes within Proteobacteria does not remain constant between seasons. The contribution of the dominant Alphaproteobacteria decreased in summer, whereas that of Betaproteobacteria increased.
Results regarding the composition of the microbial community of the Mediterranean phyllosphere present similarities with reports from other ecosystem types of the world with different assemblages of plant species. Proteobacteria is most often the dominant group, with a participation ranging from less than 40% in tropical [4] to more than 80% in temperate forests [28,29], on floating macrophytes [59] or in the phyllosphere of individual crop species, such as Spinacea oleracea [60]. Depending on the conditions prevailing, Actinobacteria, Acidobacteria, Firmicutes, Bacteriodetes are reported as second in rank phyla.
Two Sphingomonasand one Methylobacterium-related OTUs were amongst the five most abundant OTUs in the Mediterranean phyllosphere. Epiphytic microbial populations of soybean, clover and Arabidopsis [14] are also reported to be largely composed of Sphingomonas and Methylobacterium, and a Sphingomonas strain was among the five most abundant strains in the tree phyllosphere of the Brazilian Atlantic forest [27]. The dominance in the phyllosphere of methylotrophic and other taxa consuming one-carbon compounds, or of enzymes involved in the related metabolic reactions is also reported in a number of other cases [4,5,18,25,61] and has been associated with degradation processes of compounds that are toxic to plants, to humans or to the environment that are carried out by epiphytic microorganisms [5]. Dominance of these taxa in different plants, ecosystem types or geographical areas suggests major advantages from this type of symbiosis, which needs to be further explored.
It is reported that the phyllosphere microbial communities are very diverse in terms of richness, but not so in terms of evenness and that, generally, they comprise a few very well represented taxa and a large number of very rare ones [62]. This is also the case for the Mediterranean phyllosphere that we studied. Only 2% of the phyllosphere OTUs were identified as abundant, whereas 87% were identified as rare. Furthermore, the microbial community in this Mediterranean ecosystem is made up almost exclusively of habitat specialists with only very few habitat generalists.
A comparison of the epiphytic community with the airborne inoculum reveals large differences. Expressed on a per gram basis, the size of the airborne microbial community is approximately five orders of magnitude lower than that of the phyllosphere. This much thinner community bears far fewer taxa than those colonizing the phyllosphere, whereas microbes in abundance in the phyllosphere are not present in the air at the time of sampling and vice versa. A number of reasons can explain the differences between the two communities. For instance, taking air samples for a limited amount of time may not suffice to capture the full diversity of the airborne community; or, the sensitivity of the method does not allow detection at very low abundance. However, air is only a medium of microbial transfer constantly inoculating leaves with a range of taxa from various sources [63][64][65]. Microbes may land on leaves but then they are sorted out [2,66]. The epiphytic microbial community should not necessarily mirror the airborne inoculum at the time of sampling. Inocula and selection processes of the past play major roles in determining its structure.
With culture-dependent methods, it was found that the marked seasonality of the Mediterranean climate is not reflected in the size of the epiphytic microbial community [32]. Similar is the result with the culture-independent method that we used: the level of colonization did not differ between summer and winter. This was also true for Proteobacteria and Actinobacteria, which make the largest part of the community, and for several other phyla, but not for Firmicutes, Gemmatimonadetes and Chroroflexi, and the major Proteobacteria classes. In contrast, there was a seasonal effect on richness. This was higher in summer for the entire bacterial community and for most of the taxonomic groups examined, with none of the other taxa showing higher values in winter. This higher summer richness did not always lead to higher values of the diversity indices. The Pielou's evenness index was higher in summer suggesting a more homogenous OTU distribution in this season compared to winter. The Shannon diversity index differed also between seasons (higher in summer) but the Simpson index did not. This is explained by the fact that the Simpson index puts very little weight to rare species, which are numerous in the community, particularly in summer.
Experimental evidence suggests that the richness of the epiphytic bacterial community responds to climatic stressors and that it is very much influenced by drought [31]. The epiphytic microbial community of the Mediterranean phyllosphere that we studied is clearly richer in summer. This combined with the fact that the overall abundance does not change in summer, although the relative abundance of the dominant Alphaproteobacteria falls, suggest that the stressful summer conditions affect primarily the dominant members of the community. The relaxation of communities from extreme dominance of some of its members gives the opportunity to others to establish. This could be the case of Chloroflexi, which seems to be a typically summer phylum in the Mediterranean phyllosphere, although its representatives have been reported from various environments including anaerobic [67,68] and high-mountain ones [69].
The highest number of co-occurring/co-abundant OTUs that were detected through network analysis in winter corresponds mainly to Alphaproteobacteria, whereas in summer to Actinobacteria. Betaproteobacteria, with a negligible role in the winter phyllosphere network, increased considerably their co-occurrence patterns in summer. This change is in accordance with the increase in Betaproteobacteria in both richness and abundance and suggests an important role of this class in community structure in summer. Additionally, the summer network exhibits denser inter-associations among the correlated OTUs, as shown by the topological parameters calculated for this network compared to those for the winter and random networks. It has been suggested that robustness can increase with higher connectivity and denser taxa interactions, independent of taxa abundance [70], in the same way that increased taxa richness can stabilize community structure against environmental changes [71]. This is not the first time that season is found to play an important role in determining the composition of the phyllosphere bacterial community. Rastogi et al. [22] reported Proteobacteria, Firmicutes, Bacteroidetes and Actinobacteria to be the most abundantly represented phyla on lettuce foliage, as we found for the Mediterranean phyllosphere community, and a clear separation between winter and summer samples. In another study of a single Magnolia grandiflora tree, in which sampling took place at four seasons in one year and repeatedly at one season for three years [25], great temporal changes were detected, with the summer leaf community being very distinct. These temporal changes were not only seasonal; communities sampled at the same time in different years showed considerable differences, what made authors argue that seasonal patterns may not be predictable from year to year.
Knowledge of the structure of the epiphytic microbial communities in different environments and how these change with time will contribute to answering open questions on the specific functions of the microbiome on plant leaves and of the specific benefits of the partnership [72]. We found the microbial community of the Mediterranean phyllosphere to differ considerably from the air inoculum at the time of sampling, indicating selection by plants of the microbial community to be established on their leaves. This epiphytic community is dominated by habitat specialists and becomes much richer, more distinct (according to the number of OTUs on a single sample), even (according to the Pielou's eveness index) and diverse (according to the Shannon index), denser and more connected (according to network analysis) in summer. These summer features of the epiphytic microbial community could be considered as contributing to community stability under the adverse, hot and dry conditions of the Mediterranean summer, whereas the similarity in bacterial abundance in both seasons suggests that resources may not suffice for much higher population sizes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2607/7/11/518/s1, Supplement Table S1. Weather conditions prevailing on the sampling days. Measurements on the spot were taken at the end of sampling. The other values of the meteorological parameters are from the weather station at Neos Marmaras, the nearest to the sampling site, at 13 km on a straight line, which is operated by the University of the Aegean in collaboration with the Forest Office of Polygyros and the National Observatory of Athens. Supplement Table S2. OTUs (54, in total) that are identified as responsible for the within cluster similarities according to SIMPER analysis, their closest relative based on BLAST searches against SILVA 119 database, and the isolation source of the closest relative in that database. Supplement Figure S1. Rarefaction curves representing the number of OTUs against the number of high-quality reads. Funding: The research was funded through the project "ESEPMINENT", which was implemented under the "ARISTEIA II" Action of the Operational Program "Education and Lifelong Learning" and was co-funded by the European Social Fund (ESF) and National resources.

Conflicts of Interest:
The authors declare no conflict of interest.