Bacterial Communities in the Rhizosphere and Phyllosphere of Halophytes and Drought-Tolerant Plants in Mediterranean Ecosystems

The aim of the study was to investigate the bacterial community diversity and structure by means of 16S rRNA gene high-throughput amplicon sequencing, in the rhizosphere and phyllosphere of halophytes and drought-tolerant plants in Mediterranean ecosystems with different soil properties. The locations of the sampled plants included alkaline, saline-sodic soils, acidic soils, and the volcanic soils of Santorini Island, differing in soil fertility. Our results showed high bacterial richness overall with Proteobacteria and Actinobacteria dominating in terms of OTUs number and indicated that variable bacterial communities differed depending on the plant’s compartment (rhizosphere and phyllosphere), the soil properties and location of sampling. Furthermore, a shared pool of generalist bacterial taxa was detected independently of sampling location, plant species, or plant compartment. We conclude that the rhizosphere and phyllosphere of native plants in stressed Mediterranean ecosystems consist of common bacterial assemblages contributing to the survival of the plant, while at the same time the discrete soil properties and environmental pressures of each habitat drive the development of a complementary bacterial community with a distinct structure for each plant and location. We suggest that this trade-off between generalist and specialist bacterial community is tailored to benefit the symbiosis with the plant.


Introduction
The rhizosphere is considered as a microbiota hotspot that hosts a variety of organisms, comprising one of the most complicated and dynamic ecosystems on the planet [1]. Bacterial communities of the roots can be diverse and perform a number of multifaceted functions, equally important for the survival, the growth and the fitness of the plant [2], but also causing suppression of plant development under adverse conditions [3]. Indeed, rhizobacteria produce numerous compounds, which in conjunction with soil properties and environmental conditions, can trigger a dynamic system of interplay at and around the rhizosphere [4]. On the other hand, phyllospheric bacterial assemblages, being among the most ubiquitous microbial communities [5], regulate the processes in the interface between plants and Tomato rhizosphere and phyllosphere samples were collected from eight individuals in each tomato field in the Santorini Island on June 2019, grown under drought conditions. The local tomato cultivar is well-adapted to soils with volcanic properties, high light intensity and temperature, and zero irrigation scheme. No chemical pesticides for common diseases were applied for several weeks prior sampling, while basic fertilization was applied to the soil prior sowing of the seeds. Regarding the other two sampling sites, three characteristic dominant plants from each area, and five individuals from each plant were sampled. In the forest of Seich-Sou, which is characterized by acidic soils, the plants sampled were the drought-tolerant aromatic Cistus sp., Thymus sp., and Mentha pulegium (September 2018), after a period of 38 days with no rainfall in the area (climatological data of HAO DEMETER). In the National Park of Delta Axios, samples from the native halophytes Sarcoccornia sp., Atriplex sp., and Crithmum sp. were collected in June 2018. This wetland covers >33,000 hectares of land, includes estuaries of four rivers and has been included in the Natura 2000 network of European ecological regions [29]. The sampling site in the wetland was located proximate to the estuary of the River Axios, an either saline or sodic environment, characterized by alkaline soils. All plant samples were placed in sterile bags and brought back to the lab under sterile and cold conditions within 6 h.

Soil Properties
Soil samples were sampled from five different points in each site, and mixed in the same ratio to form a compiled sample. Soil samples were then air-dried, ground to pass a 2-mm sieve and analysed for particle size distribution according to Bouyoucos [30], and for chemical properties according to Sparks [31]. Briefly, pH was measured in a (1:2) soil suspension with water. The electrical conductivity (EC se ) and the water-soluble cations K, Na, Ca and Mg were determined in the saturation extract. For the analytical determinations of K and Na, flame photometry was used, whereas Ca and Mg were Microorganisms 2020, 8, 1708 4 of 20 determined by atomic absorption spectrometry. Sodium adsorption ratio (SAR) was calculated from the concentrations of the water-soluble Na, Ca and Mg. All the above methods have been described analytically elsewhere [32]. Furthermore, organic C was determined by the wet oxidation method and CaCO 3 with a volumetric calcimeter [33]. Soil available NO 3 -N was extracted with 1 M KCl and determined by ultraviolet spectrometry [34]. Olsen-P was extracted with 0.5 M NaHCO 3 , pH 8.5 and determined by the molybdenum blue-ascorbic acid method [35]. Exchangeable K was extracted with 1 M CH 3 COONH 4 , pH 7 and determined by flame photometry [36].

Sample Processing and Sequencing
To isolate the bacterial communities of the rhizospheric and phyllospheric samples, the bulk soil, and other external material was removed with manual shaking of the root and leaves, and 0.5-1 g from each sample was transferred into phosphate saline buffer (PBS; NaCl 137 nmol L −1 , KH 2 PO 4 1.8 nmol L −1 , KCl 2.7 nmol L −1 and Na 2 HPO 4 1.42 nmol L −1 , pH = 7.4) and sonicated for 10 min (Transsonic 460). The solutions were subsequently centrifuged at 10000 rpm for 20 min, and the sedimentation material was placed in −80 • C until DNA extraction. The rhizosphere samples consisted of both the primary and lateral roots. DNA was extracted from a total of 92 samples (46 rhizospheric and 46 phyllospheric from all plants in all areas), using a NucleoSpin ® Soil Genomic DNA Isolation Kit (Macherey-Nagel, Bethlehem, PA, USA), according to the manufacturer's instructions. The concentration and quality of the recovered DNA was checked by a NanoDrop TM spectrophotometer (Thermo Scientific TM , Waltham, MA, USA) and confirmed via gel electrophoresis.
The extracted DNA was subjected to PCR using the specific primers targeting the V6-V8 hyper-variable region of the 16S rRNA gene (B969F = ACGCGHNRAACCTTACC and BA1406R = ACGGGCRGTGWGTRCAA) [37]. These primers have been found to successfully amplify approximately 470 bp of all the major high-level bacterial taxonomic groups with up to 83% coverage, after in silico analysis via the SILVA TestPrime 1.0, performed at the Integrated Microbiome Resource (IMR) at Dalhousie University (Halifax, NS, Canada). The PCR products were verified visually by running on a high-throughput Hamilton Nimbus Select robot (Hamilton Company, Reno, NV, USA) using Coastal Genomics Analytical Gels (Hamilton Company) and were cleaned-up and normalized using the high-throughput Charm Biotech Just-a-Plate 96-well Normalization Kit (Charm Biotech, Lawrence, MA, USA). The amplicon samples were sequenced on Illumina MiSeq using 300 + 300 bp paired-end chemistry which allows for overlap and stitching together of paired amplicon reads into one full-length read (http://cgeb-imr.ca/protocols.html). The PCR amplification step and sequencing steps were performed at the Integrated Microbiome Resource (IMR) at Dalhousie University (Halifax, NS, Canada).

Data Analysis
All analyses were performed on the rarefied dataset. Alpha-diversity estimators (the richness estimator S Chao1 , the Simpson, and Equitability indexes) were calculated for all samples with the PAST 3.16 software [38]. One-way analysis of similarities (ANOSIM) based on Bray-Curtis dissimilarity coefficient was used to examine whether the similarity in bacterial composition of the phyllosphere and rhizosphere between sites was greater than or equal to the similarity within the sites. Then, we used one-way permutational multivariate analysis of variance (PERMANOVA) based on the Bray-Curtis dissimilarity coefficient in the PAST 3.16 software [38] to test whether the soil properties separate the rhizospheric bacterial assemblages of the different sampling sites. To compute beta diversity of the plants' bacterial assemblages for each sampling site for the rhizosphere and phyllosphere separately, we used the 'betapart' R package version 1.5.1 [39]. Beta diversity is a measure of the difference in species composition between two or more local assemblages [40], partitioned into two components according to Baselga et al. [39,41]: (i) spatial turnover in species composition, measured as Simpson dissimilarity (bSIM); and (ii) variation in species composition due to nestedness (bNES) measured as nestedness-resultant fraction of Sorensen dissimilarity. Based on the above, there are two potential Microorganisms 2020, 8, 1708 5 of 20 ways in which two or more assemblages can be different: one is species replacement (turnover) and the second is species loss or gain, which implies that the poorest assemblage is a subset of the richest one (nestedness). Then, the bacterial assemblages of the different samples (both rhizosphere and phyllosphere samples) were compared using the Plymouth routines in the multivariate ecological research software package PRIMER v.6 [42]. The Jaccard similarity coefficients were calculated based on OTUs presence/absence data, to identify interrelationships between samples and construct an nMDS plot.
Furthermore, OTUs were classified as abundant when their total number of reads exceeded the 0.2% of the total number of reads of the entire dataset (i.e., >900 reads in all samples in these OTUs). In order to identify OTUs with ubiquitous presence in all plants and locations, pointing out to a common bacterial community shared between plants regardless of soil properties, plant characteristics, plant compartment examined, location and environmental conditions, the Levins' index was calculated [43]. 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. [44], using Levins' niche width (B) index [43]: where Pij is the proportion of OTU j in sample i, and N is the total number of samples. Therefore, B describes the extent of niche specialization based on the distribution of OTU abundances without considering 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 33 (top generalist). The OTUs with B index higher than 20 were arbitrarily considered as generalists, while the OTUs with B lower than 10 were categorized as specialists [45].

Read Processing
The produced reads were subjected to downstream processing using the mothur v.1.34.0 software [46], following the proposed standard operating procedure [47]. Briefly, forward and reverse reads were joined, and contigs below 200 bp, with >8 bp homopolymers and ambiguous base calls were removed from further analysis. The remaining reads were dereplicated to the unique sequences and aligned independently against the SILVA 132 database, containing 1,861,569 bacterial SSU rRNA gene sequences [48]. The reads suspected of being chimeras were removed using the UCHIME software [49]. The remaining reads were clustered into Operational Taxonomic Units (OTUs) at 97% similarity level, using the average neighbor method in mothur. 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 [50,51]. The resulting dataset was rarefied with the subsample command in mothur v.1.34.0 to 7125 reads, while 11 samples with lower total number of reads were also included in the dataset. We chose this course of action as the best compromise in order to retrieve the majority of the biodiversity detected by sequencing, as rigorous subsampling may result in some OTUs loss, but also to include all samples in the analysis, even those with lower number of reads. Overall, 77 samples (of the 92) were retained for further data analysis, while 15 were excluded because of low sequencing success (Table 1). Taxonomic classification was assigned using BLASTN searches against the SILVA 132 [48], with curated bacterial taxonomy, by applying the lowest accepted level of >80% similarity with a closest relative. The reads belonging to OTUs related to Chloroplasts were removed from the dataset. The raw reads were submitted to GenBank-SRA under the BioProject accession number PRJNA635261.

Soil Properties in the Sampling Sites
Detailed description of the soil properties in the sampling sites is reported elsewhere [52]. Briefly, the soils collected from the tomato cultivar in the two sampling sites of Santorini (Vlichada and Emporio) were inorganic, coarse in texture, calcareous and alkaline in reaction. They represent ecosystems with different soil fertility based on the contained available NO 3 -N, Olsen P and K: In Emporio, these concentrations were similar or higher than the critical sufficiency range, i.e., 10 mg kg −1 NO 3 -N, 10 mg kg −1 P and 110 mg kg −1 K, respectively, whereas the opposite was evident for Vlichada, in which the concentrations of the three macronutrients were almost five times lower than those of Emporio. Even though they were considered as nutrient-poor and nutrient-rich sites, respectively, they were clustered as more similar compared to the other sampling sites according to Euclidean distances ( Figure 1). The soil in the Seich-Sou forest was acidic, enriched with organic matter, and contained adequate amounts of available NO 3 -N and P, but low concentrations of K + . On the contrary, the soil from the National Park of Delta Axios was completely different according to Euclidean distances of its properties (Figure 1), and it was an alkaline (pH = 8.1) saline-sodic soil. This was evidenced by the values of Electrical Conductivity (EC) = 69 mS cm −1 and SAR 81 [52]. Overall, Vlichada and Emporio were clustered together in terms of Euclidean distances of their soil properties, and the other two locations were dissimilar, with closer to the Santorini sampling sites being the Seich-Sou forest ( Figure 1). to Chloroplasts were removed from the dataset. The raw reads were submitted to GenBank-SRA under the BioProject accession number PRJNA635261.

Soil Properties in the Sampling Sites
Detailed description of the soil properties in the sampling sites is reported elsewhere [52]. Briefly, the soils collected from the tomato cultivar in the two sampling sites of Santorini (Vlichada and Emporio) were inorganic, coarse in texture, calcareous and alkaline in reaction. They represent ecosystems with different soil fertility based on the contained available NO3-N, Olsen P and K: In Emporio, these concentrations were similar or higher than the critical sufficiency range, i.e., 10 mg kg −1 NO3-N, 10 mg kg −1 P and 110 mg kg −1 K, respectively, whereas the opposite was evident for Vlichada, in which the concentrations of the three macronutrients were almost five times lower than those of Emporio. Even though they were considered as nutrient-poor and nutrient-rich sites, respectively, they were clustered as more similar compared to the other sampling sites according to Euclidean distances ( Figure 1). The soil in the Seich-Sou forest was acidic, enriched with organic matter, and contained adequate amounts of available NO3-N and P, but low concentrations of K + . On the contrary, the soil from the National Park of Delta Axios was completely different according to Euclidean distances of its properties (Figure 1), and it was an alkaline (pH = 8.1) saline-sodic soil. This was evidenced by the values of Electrical Conductivity (EC) = 69 mS cm −1 and SAR 81 [52]. Overall, Vlichada and Emporio were clustered together in terms of Euclidean distances of their soil properties, and the other two locations were dissimilar, with closer to the Santorini sampling sites being the Seich-Sou forest ( Figure 1).

Bacterial Communities' Diversity and Composition
The sequencing pipeline was proven more effective in the rhizospheric samples; it is noteworthy that only two rhizospheric samples were excluded after data processing because of low overall number of sample reads (from Cistus sp. and Thymus sp. samples in the Seich-Sou forest), while 13 phyllospheric samples (from the Seich-Sou forest and the National Park of Delta Axios) had low sequencing success and thus removed from the dataset (Table 1). For further analysis, 44 out of 46 samples obtained from the rhizosphere (96% of the total rhizosphere samples) and 33 out of the 46 samples from the phyllosphere (72% of the samples) were kept. Rarefaction curves were constructed for all samples, indicating a good coverage of the bacterial communities' diversity in most of them, and in at least one sample from each plant (Supplementary Figure S2). The ratio of observed to expected (S Chao1 ) number of OTUs in all samples was >0.55, but in at least one sample per plant the ratio reached >0.75 (Supplementary Table S1

Bacterial Communities' Diversity and Composition
The sequencing pipeline was proven more effective in the rhizospheric samples; it is noteworthy that only two rhizospheric samples were excluded after data processing because of low overall number of sample reads (from Cistus sp. and Thymus sp. samples in the Seich-Sou forest), while 13 phyllospheric samples (from the Seich-Sou forest and the National Park of Delta Axios) had low sequencing success and thus removed from the dataset (Table 1). For further analysis, 44 out of 46 samples obtained from the rhizosphere (96% of the total rhizosphere samples) and 33 out of the 46 samples from the phyllosphere (72% of the samples) were kept. Rarefaction curves were constructed for all samples, indicating a good coverage of the bacterial communities' diversity in most of them, and in at least one sample from each plant (Supplementary Figure S2). The ratio of observed to expected (SChao1) number of OTUs in all samples was >0.55, but in at least one sample per plant the ratio reached >0.75 (Supplementary Table S1 Table 1 for coding assignments. Overall, 80% of the total number of OTUs recovered from all samples were affiliated to 21 highlevel taxonomic groups and the rest were affiliated to unclassified/unidentified Bacteria (data not shown). Among the 21 groups that were detected, eight of them comprised of >70% of the total number of OTUs ( Figure 3). In all samples, Proteobacteria was the dominant group, comprising between 33 and 47% of the total OTUs richness in each plant (either in the rhizosphere or in the phyllosphere), followed by Actinobacteria (8-13% of the total number of OTUs in all samples) and Bacteroidetes (6-11%). Within Proteobacteria, Alphaproteobacteria was the group with the highest OTUs richness (16-22%), followed by Gammaproteobacteria (6-10%) and Deltaproteobacteria (5-7%) (Figure 3). Other less diverse taxonomic groups included the groups Chloroflexi, Deinococcus, Nitrospirae, Deinococcus-Thermus and others with <10 OTUs in the dataset.  Table 1 for coding assignments.
Overall, 80% of the total number of OTUs recovered from all samples were affiliated to 21 high-level taxonomic groups and the rest were affiliated to unclassified/unidentified Bacteria (data not shown). Among the 21 groups that were detected, eight of them comprised of >70 % of the total number of OTUs (Figure 3). In all samples, Proteobacteria was the dominant group, comprising between 33 and 47% of the total OTUs richness in each plant (either in the rhizosphere or in the phyllosphere), followed by Actinobacteria (8-13% of the total number of OTUs in all samples) and Bacteroidetes (6-11%). Within Proteobacteria, Alphaproteobacteria was the group with the highest OTUs richness (16-22%), followed by Gammaproteobacteria (6-10%) and Deltaproteobacteria (5-7%) (Figure 3). Other less diverse taxonomic groups included the groups Chloroflexi, Deinococcus, Nitrospirae, Deinococcus-Thermus and others with <10 OTUs in the dataset.  Table 1 for coding assignments.

Rhizosphere vs. Phyllosphere and Spatial Heterogeneity
Rhizospheric samples were found to be more diverse than phyllospheric ones in all locations and plants, with higher overall OTUs richness in all cases. In addition, commonly found OTUs between the rhizosphere and the phyllosphere of the same plant represented only a portion of the total number of OTUs detected in either plant compartment; overall only 3735 out of the 16,355 were common between the two plant compartments (Figure 4). Evidently, in most of the cases the phyllospheric samples per plant that were included in the analysis were fewer than the rhizospheric ones. In the case of the tomato plants of Vlichada and Emporio and Thymus sp. of the Seich-Sou forest (which included the same number of analysed samples for each compartment), the number of OTUs detected in the rhizosphere was much higher (Figure 4). Taking into consideration that the rhizospheric and phyllospheric samples appeared to consist of heterogenous bacterial assemblages between them, the beta diversity of the bacterial communities of all plants in one location was for  Table 1 for coding assignments.

Rhizosphere vs. Phyllosphere and Spatial Heterogeneity
Rhizospheric samples were found to be more diverse than phyllospheric ones in all locations and plants, with higher overall OTUs richness in all cases. In addition, commonly found OTUs between the rhizosphere and the phyllosphere of the same plant represented only a portion of the total number of OTUs detected in either plant compartment; overall only 3735 out of the 16,355 were common between the two plant compartments (Figure 4). Evidently, in most of the cases the phyllospheric samples per plant that were included in the analysis were fewer than the rhizospheric ones. In the case of the tomato plants of Vlichada and Emporio and Thymus sp. of the Seich-Sou forest (which included the same number of analysed samples for each compartment), the number of OTUs detected in the rhizosphere was much higher (Figure 4). Taking into consideration that the rhizospheric and phyllospheric samples appeared to consist of heterogenous bacterial assemblages between them, the beta diversity of the bacterial communities of all plants in one location was for both compartments higher than 0.85 regardless the sampling site, attributed to turnover of >0.77. The nestedness remained low with the lowest values being computed for the plants in the National Park of Delta Axios (Table 2).
phyllosphere bacterial communities. According to one-way ANOSIM the similarity of the bacterial composition within the same site was greater for both the phyllosphere (R = 0.58, psame < 0.001) and the rhizosphere (R = 0.91, psame < 0.001) than between the different sites. The significant differences between the site pairs for both the phyllosphere and the rhizosphere was supported by one-way PERMANOVA which detected significant effects of soil properties (F = 4.43, p < 0.001 for the phyllosphere; and F = 12.9, p < 0.001 for the rhizosphere) on the bacterial community structure.  Table 1 for coding assignments. Furthermore, the Jaccard similarity index showed a clear separation of the bacterial assemblages according to the sampling site ( Figure 5), thus corroborating PERMANOVA indications on soil properties/location effects on bacterial community structure. Apart from the bacterial assemblages of the tomato plants in Emporio, Santorini, the rhizosphere samples were clearly separated from the phyllosphere samples, further confirming the results of beta diversity estimates and the one-way ANOSIM, that rhizosphere bacterial communities are dissimilar to the respective phyllosphere  Table 1 for coding assignments. According to one-way ANOSIM the similarity of the bacterial composition within the same site was greater for both the phyllosphere (R = 0.58, p same < 0.001) and the rhizosphere (R = 0.91, p same < 0.001) than between the different sites. The significant differences between the site pairs for both the phyllosphere and the rhizosphere was supported by one-way PERMANOVA which detected significant effects of soil properties (F = 4.43, p < 0.001 for the phyllosphere; and F = 12.9, p < 0.001 for the rhizosphere) on the bacterial community structure.
Furthermore, the Jaccard similarity index showed a clear separation of the bacterial assemblages according to the sampling site ( Figure 5), thus corroborating PERMANOVA indications on soil properties/location effects on bacterial community structure. Apart from the bacterial assemblages of the tomato plants in Emporio, Santorini, the rhizosphere samples were clearly separated from the phyllosphere samples, further confirming the results of beta diversity estimates and the one-way ANOSIM, that rhizosphere bacterial communities are dissimilar to the respective phyllosphere communities. Finally, smaller groupings of the samples belonging to the same plant/same compartment were also apparent ( Figure 5). communities. Finally, smaller groupings of the samples belonging to the same plant/same compartment were also apparent ( Figure 5).  Table 1 for coding assignments.

Generalists and Abundant Specialists
According to Levins' index (B), 36 OTUs were found to be generalists (Table 3), either by taking into consideration the entire dataset as a whole, or by calculating B in each location, in the rhizosphere and phyllosphere datasets separately and identifying the common generalists in all cases. Euclidean distances based on the number of reads in each generalist OTU in all samples, showed similar groupings as indicated by the Jaccard similarity index (Figure 6), separating sampling sites (thus soil properties), plant species and rhizospheric vs. phyllospheric samples. Most of the generalists belonged to Actinobacteria (12 OTUs out of the 36), followed by Alphaproteobacteria (8/36). Generalist OTUs were also assigned to Bacteroidetes (1 OTU of the 36), Firmicutes (1/36) and Gammaproteobacteria (1/36), while 13 OTUs were assigned to an uncultured bacterial strain (Table 3).  Table 1 for coding assignments.

Generalists and Abundant Specialists
According to Levins' index (B), 36 OTUs were found to be generalists (Table 3), either by taking into consideration the entire dataset as a whole, or by calculating B in each location, in the rhizosphere and phyllosphere datasets separately and identifying the common generalists in all cases. Euclidean distances based on the number of reads in each generalist OTU in all samples, showed similar groupings as indicated by the Jaccard similarity index (Figure 6), separating sampling sites (thus soil properties), plant species and rhizospheric vs. phyllospheric samples. Most of the generalists belonged to Actinobacteria (12 OTUs out of the 36), followed by Alphaproteobacteria (8/36). Generalist OTUs were also assigned to Bacteroidetes (1 OTU of the 36), Firmicutes (1/36) and Gammaproteobacteria (1/36), while 13 OTUs were assigned to an uncultured bacterial strain (Table 3).
On the other hand, 22 OTUs were characterized as abundant specialists, meaning they appeared in high abundances in individual samples and thus appeared to be overall abundant with >900 number of reads in the entire dataset, while they were rare or absent in the majority of the samples (Table 4). Among these, 10 belonged to Alphaproteobacteria, six to Gammaproteobacteria, four to Bacteroidetes, one to Firmicutes and one to Actinobacteria. The specialist nature of these abundant OTUs was evident by the habitats of preference, where they increased their number of reads by >20% of the sample's reads. In the majority of cases they were detected in only one habitat: e.g., OTU_2 and OTU_5/OTU_66, in Vlichada and Emporio site of Santorini Island, respectively; OTU_10, OTU_13, OTU_15, OTU_28, OTU_38, OTU_93, OTU_133, OTU_142 and OTU_144 in the Seih-Sou forest; OTU_7, OTU_9, OTU_22, OTU_37, OTU_61, OTU_68, OTU_103 and OTU_149 in the National Park of Delta Axios. By contrast, only in few cases abundant specialists were detected in certain plant species (e.g., OTU_2, OTU_5 and OTU_66 in tomatoes of Santorini; OTU_10 in Mentha pulegium, OTU_15, OTU_28, OTU_142 and OTU_144 in Cistus sp., OTU_5 in Thymus sp. of the Seich-Sou forest; OTU_22, OTU_37, OTU_103 in Sarcocornia sp., OTU_68 and OTU_149 in Crithmum sp. of the National Park of Delta Axios) ( Table 4). Table 3. The OTUs identified as being rhizosphere generalists in the three sampling areas (Santorini Island, Seich-Sou Forest and the National Park of Delta Axios), their putative high-level taxonomic affiliation, their closest relatives based on BLAST searches against the SILVA database and confirmed against NCBI database, and the isolation source of the closest relative. communities. Finally, smaller groupings of the samples belonging to the same plant/same compartment were also apparent ( Figure 5).  Table 1 for coding assignments.

Generalists and Abundant Specialists
According to Levins' index (B), 36 OTUs were found to be generalists (Table 3), either by taking into consideration the entire dataset as a whole, or by calculating B in each location, in the rhizosphere and phyllosphere datasets separately and identifying the common generalists in all cases. Euclidean distances based on the number of reads in each generalist OTU in all samples, showed similar groupings as indicated by the Jaccard similarity index (Figure 6), separating sampling sites (thus soil properties), plant species and rhizospheric vs. phyllospheric samples. Most of the generalists belonged to Actinobacteria (12 OTUs out of the 36), followed by Alphaproteobacteria (8/36). Generalist OTUs were also assigned to Bacteroidetes (1 OTU of the 36), Firmicutes (1/36) and Gammaproteobacteria (1/36), while 13 OTUs were assigned to an uncultured bacterial strain (Table 3).  Table 4. The OTUs identified as being abundant specialists in the entire dataset, their putative high-level taxonomic affiliation, their closest relatives based on BLAST searches against the SILVA database and confirmed against NCBI database, and the isolation source of the closest relative. Shaded boxes denote high relative abundance; i.e >20% of the sample's total number of reads on average of total sampled individuals in each occasion and different shading colours represent different locations; i.e., black denotes Santorini Island; grey, Seih-Sou forest; light grey, the National Park of Delta Axios.

Discussion
In the present study, the diversity and variability of bacterial communities in the rhizosphere and phyllosphere of native wild plants and local tomato cultivars of three different diversely stressed Mediterranean ecosystems were examined via 16S rRNA gene amplicon sequencing. Overall, in all sampling locations, a large bacterial diversity was revealed. In almost all plant species, the rhizospheric bacterial diversity, in terms of species richness, was much higher than the respective phyllopsheric communities. This is a common finding in similar studies of native and/or cultivated plants in different environments [53][54][55][56][57][58][59][60]. These differences in species richness between the two plant compartments have been attributed to their fundamental physiological and functional differences, as well as the direct impact of their surrounding environments (soil vs air environment). It is well established that the root exudates contribute to the selective growth of specific bacteria [19,56] by signalling plant-microbe and microbe-microbe interactions [61], and can promote the differentiation of the bacterial assemblages through soil-driven selection [59,60]. On the other hand, phyllosphere, being a harsh habitat with shorter lifespan, less nutrient availability [5], and characterized by swift fluctuations in environmental pressures [62], exhibits different physiological functions [18] leading to generally lower bacterial richness and abundance. In the present study, a marked higher richness of bacterial OTUs in the rhizosphere was apparent, reflecting the influence of the stressed habitats and the variability of mechanisms developed by the sampled plants and other co-existing organisms to exploit all available synergistic survival features of the surrounding soil environment. The fact that high-quality DNA was extracted and analysed in more rhizospheric samples, possibly contributing to the higher rhizobacterial OTUs richness, might also be an additional indication of low phyllosphere bacterial abundance.
Our results further indicated that the bacterial communities varied significantly in the three sampling locations. In particular, bacterial communities of the plants in the National Park of Delta Axios, which consists of alkaline and saline-sodic soils, were grouped separately from the samples of the other two locations, according to the Jaccard index, regardless of the plant species and/or the plant compartment. In contrast, bacterial assemblages from the Seich-Sou forest and the Santorini Island exhibited similarities in respect to sampling site (10% similarity, Jaccard index), and stronger on species level within the same site (20% similarity, Jaccard index), probably reflecting similar same soil properties and environmental characteristics. These results are in accordance with the general notion that soil microbial biogeography is primarily controlled by edaphic properties [63], which consequently affect the rhizobacterial communities' structure in different soils. The bacterial communities of the tomato cultivars of the Santorini Island, on both locations of Vlichada and Emporio, characterized as relatively nutrient-poor and nutrient-rich respectively, showed similarities in the rhizosphere. Local and regional soil properties are identified factors which strongly govern the bacterial community structure of rhizobacteria [64], as they formulate specialized niches favouring the growth of unique bacterial assemblages depending on these properties [63,65]. Indeed, soil properties of the Santorini Island are influenced by the volcanic environment of the island [66], possibly showing convergence and similar soil properties throughout the island, supporting the data of similarities in the rhizospheric communities of the two tomato cultivar fields, independently of the level of soil fertility. However, the bacterial community of the phyllosphere in Vlichada formed a distinct group, suggesting that phyllospheric bacterial assemblages are affected not only by the site, atmospheric conditions and possible air-dispersals [18,67], but also by other factors, such as nutrient availability. As Mello·et al. [68] argued, nutrient limitations may have a significant selective pressure on the biodiversity of the microorganisms present in a harsh environment. Indeed, in the nutrient-poor site of Vlichada, shortage of nutrients seemed to affect the phyllospheric communities more than the plant compartment and/or specific location, based on the Jaccard index. These findings need further investigation, taking into consideration that microbe-microbe interactions and within plant compartments microbes' dispersal are important selective forces forming the microbiomes' structure of the rhizosphere, phyllosphere and plant endosphere compartments [22].
Although overall the bacterial communities' composition and structure were similar (in terms of high taxonomic group composition) in all locations, plant species, and plant compartments, dominated mainly by Proteobacteria, followed by Actinobacteria, the 36 generalist OTUs as identified by Levins' index, belonged mostly to Actinobacteria (>30%). Several studies have shown the dominance of Proteobacteria and Actinobacteria in rhizospheric bacterial communities [69,70]. These groups represent ubiquitous rhizospheric taxa detected in various stressed environments, with many biotechnological applications in sustainable agriculture [71]. Proteobacteria and Actinobacteria were among the top generalists in the phyllospheres of nine perennial plants in a Mediterranean ecosystem [72]. Among the generalist OTUs, taxa closely related to the genera Bradyrhizobium, Steroidobacter, Arthrobacter, Mycobacterium, Pseudonocardia, Rhizobium, Bosea, and Paenibacillus have been associated with higher nutrient uptake, antioxidant activity [73], catalase activity [74], plant growth [75], phosphate solubilization [76], siderophore production [77], nitrogen fixation [78], cytokinin production [79], and other PGP traits. This consistency points towards a common core of a bacterial community present in different stressed environments of the Mediterranean region, which may establish beneficial interactions with the host plants in a variety and complementary ways.
On the other hand, OTUs that were characterized as abundant, but were found with high number of reads only in few individual samples and were rare or absent from the majority of the samples, were characterized as abundant specialist taxa and comprised of a different taxonomic composition than generalist OTUs. It has been suggested that specialization traits can successfully differentiate bacterial communities [45]. While generalists are able to adapt to varying environmental conditions and thrive in multiple habitats constituting the core of bacterial communities, specialist taxa decrease dramatically in abundance, or even disappear following minor environmental changes [80], thus leading to species sorting, filtered by local environmental conditions [81]. Most abundant specialists of the present study have been detected in high abundances (>20% of the total sample reads) in a single habitat. For example, OTUs related to the taxa Hymenobacter, Methylobacterium, Novosphingobium, and Phenylobacterium have been found in high abundances in the drought-stressed system of the Seich-Sou forest. These taxa can also be characterized as indicator taxa, when found in high abundances, for ecosystems with similar soil properties as the Seich-Sou forest [82]. Hymenobacter species are UV-resistant bacteria isolated from drought-stressed areas with high UV radiation and low temperatures [83], endosymbiotic Methylobacterium mitigates the impact of limited water availability in crops [84], Novosphingobium spp. have been associated with dry systems [85], and Phenylobacterium spp. have been reported to be favoured by heat [86]. On the other hand, the taxa Marinobacter, Gracilimonas, Lewinella and Jannaschia have been detected in >20% of the sample's reads only in the saline environment of the National Park of Delta Axios, and similarly can be characterized as indicator taxa for similar ecosystems, when in high abundances. These taxa have all been commonly isolated from marine sediments and high-salinity environments [87][88][89][90] and are considered characteristic halotolerant bacteria. Stressed environments overall can play critical roles in structuring the plant bacterial communities by selecting stress tolerant groups of microbes not necessarily correlated directly to the vegetation of the ecosystem [82]. Moreover, it seems that in the studied ecosystems, severity of unfavourable conditions and nutrient availability are both equally important aspects that shape the bacterial communities in terms of richness, structure and behavioural characteristics, which in turn affect plant-bacteria co-existence.

Conclusions
Knowledge of the composition and structure of the rhizo-and phyllospheric bacterial communities in plants of stressed environments will shed some further light on the richness, variability, dispersal, establishment, functioning of the microbiome on the plant compartments and the specific benefits of this partnership. Our results of 16S rRNA gene amplicon sequencing from rhizo-and phyllospheric samples of native plants in high-saline, dry Mediterranean ecosystems, suggested the presence of overall highly diverse, variable bacterial communities, which differed depending on the plant's compartment, the soil properties and location of sampling, pointing towards abiotic filtering and environmental drivers. In addition, a commonly found pool of generalist taxa was detected independently of sampling location, plant, or plant compartment. We suggest that these shared OTUs contribute to the plants' survival and growth through mechanisms for a wide-range and generic management of environmental stresses and complement the variable communities in each soil type which address specific stresses caused because of the plant's location. The investigation of the driving forces that shape the composition and structure of the rhizo-and phyllospheric bacterial communities in plants under stress will shed further light on the richness, variability, dispersal, establishment, and functioning of the plant microbiome, the role of generalist and specialist taxa, and the specific benefits of plant-bacteria symbioses under harsh conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2607/8/11/1708/s1, Figure S1: Sampling locations of tomato plants from two different locations in Santorini Island (Vlichada and Emporio), of aromatic plants from the forest of Seich-Sou and of halophytes from the National Park of Delta Axios. Figure S2: Rarefaction curves representing the number of OTUs against the number of high-quality reads after rarefication and removal of ambiguous reads. Table S1: Number of OTUs, the richness estimator (S Chao1 ), the ratio observed/expected OTUs and the heterogeneity of the alpha-diversity indexes (the Simpson dominance and Equitability indexes) and the total number of reads in each sample. The total and average number of OTUs, number of reads and average values of the alpha diversity estimators are also shown. N/A: Not Applicable.