Sampling was performed in April in the year 2016. Plant samples were collected at five distinct successional stages in the salt marsh chronosequence. This chronosequence spans more than 100 years of primary succession and is located on the island of Schiermonnikoog (53°30′N, 6°10′E), Netherlands (for detailed information on sampling, see Wang et al., 2016). Samples were collected at locations with successional ages of 5, 15, 35, 65, and 105 years. In this study, we use soil stage as a proxy of soil type, as previous work has revealed that the sedimentation caused by the tidal regime has resulted in modifications on the soil physicochemical conditions along the primary succession, leading to an accumulation of silt and clay particles. In addition, the salinity level also increased over time during succession, due to an accumulative effect. For the details on the establishment of sampling plots and descriptions on the chronosequence verification, see Olff et al. (1997), Dini-Andreote et al. (2014), and Wang et al. (2016). Briefly, triplicated plots within each of the locations were established at the same base of elevation [vertical position relative to mean sea level at the initial elevation gradient on the bare sand flats with a base elevation of 1.16 m ± 2.2 cm (mean ± SE) above Dutch Ordinance Level]. Within each plot, four healthy-looking L. vulgare and A. maritima of similar sizes with attached soil adhering to the intact roots were collected and processed together generating two composite samples per plot. Thirty composite samples in total were collected (5 stages × 3 plots per stage × 2 plant species). Each sample was placed into a sterile plastic bag, sealed and transported to the laboratory in <24 h. From each composite sample, we separated rhizosphere soil from plant roots (see below).

A detailed description for the workflow of the isolation and screening of plant-associated bacteria is given in Supplementary Figure S1A. Rhizosphere soil samples were collected by weighting ten grams of roots with tightly adhering soil particles (about 3 g rhizosphere soil). Root samples were transferred into an Erlenmeyer flask containing 47 mL of sterile 1X Phosphate Buffered Saline (PBS buffer) and shaken for 30 min at 200 rpm at room temperature. An aliquot of 1 mL of the suspension with rhizosphere soil was transferred into sterile 1X PBS buffer and serial dilutions (1/10) were prepared.

Plant roots (about 8 g) were thoroughly washed with running tap water, trimmed to remove adhering soil and dead tissues, followed by surface sterilization, which consisted of immersion in 1.5% NaClO solution for 3 min, followed by 70% ethanol for 3 min and sterile distilled water (3 × 3 min). The surface-sterilized root samples (5 g) were diced with a sterile scalpel and immersed into 45 mL of 0.9% NaCl solution. After shaking incubation for 1 h at 28°C, the suspension with root pieces was shaken using a horizontal vortex instrument (4 × 1 min, 30 s in-between). An aliquot of 1 mL of the suspension containing the released plant endophytes was transferred into sterile 1X PBS buffer and serial dilutions (1/10) were prepared. Sterility checks were performed by tissue-blotting surface-sterilized root samples on R2A plates at 28°C for 2–7 days. Only samples without bacterial growth were considered successfully sterilized and used further.

R2A medium, as a non-selective medium recommended for the examination of total heterotrophic bacteria in soil (Ellis et al., 2003), was used to culture the heterotrophic population. Aliquots of 0.1 mL of the 10^-1 to 10^-3 dilutions from rhizosphere soil and root samples were spread on R2A medium plates and incubated for 2 days at 25°C, after which we determined the number of colony forming units (CFU). A maximum of 32 bacterial colonies with unique morphologies per replicate were purified using a streak-plate procedure (see Supplementary Table S1 for total numbers), transferred onto new R2A medium plates and further used as templates for BOX-PCR – a DNA-based typing method that differentiate bacterial species at strain level by simultaneously screening DNA regions scattered in the bacterial genome (Brusetti et al., 2008). To improve cell lysis, the colonies were first inoculated into 50 μL of NaOH solution (0.05 M) and then lysed at 95°C for 15 min in the PCR machine. BOX-PCR was performed by using the BOX-A1R primer (5′-CTACGGCAAGGCGACGCTGACG-3′) (Versalovic et al., 1994). Twenty μL PCR reactions were performed using 0.32 μL 5 U μL^-1 Taq DNA Polymerase, 4 μL of 5X Gitschier Buffer [83 mM (NH₄)₂SO₄, 335 mM Tris-HCl (pH 8.8), 32.5 mM MgCl₂, 325 mM EDTA (pH 8.8), 1% commercial stock of β-mercaptoethanol, ddH2O], 2 μL 100% DMSO, 1 μL 25 mM of each dNTP in a mixture, 0.32 μL 20 mg mL^-1 bovine serum albumin (BSA) (Roche Diagnostics GmbH, Mannheim, Germany), and 0.8 μL of 10 μM BOX-A1R primer and 1 μL of the lysed cell solution. The thermal cycler protocol was 95°C for 3 min, 35 cycles of 94°C for 4 s, 92°C for 30 s, 50°C for 60 s, 65°C for 8 min and a final 16 min extension at 65°C. BOX-PCR profiles were visualized by separation on 2% agarose gel and staining with ethidium bromide. Images of the gels were visualized and documented under UV light with Image Master VDS system (Amersham Biosciences, United Kingdom). For the details of the bacterial isolates showing unique BOX-PCR profiles, see Supplementary Tables S1 and S2. Culture stocks for individual isolate were stored in 25% glycerol at -80°C.

A total of 159 bacterial cultures with unique BOX-PCR patterns (Supplementary Figure S1 and Table S2) were subjected to total DNA extraction using the MoBio UltraClean Microbial DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, United States). We followed the instruction manual, except for heating the preparations at 65°C for 10 min with occasional bump vortexing for a few seconds every 2–3 min. The amount of DNA in each sample was quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies). All DNA samples were standardized to the equal concentration of 5 ng μL^-1 for further analyses.

16S rRNA gene amplification was performed by using the bacterial-specific primers, B8F (5′-AGAGTTTGATCMTGGCTCAG-3′) (Edwards et al., 1989) and U1406R (5′-ACGGGCGGTGTGTRC-3′) (Lane, 1991). Fifty μL PCR reactions were performed with 0.2 μL 5U μL^-1 FastStart High Fidelity (FSHF) Taq DNA Polymerase, 5 μL 10X FSHF Reaction buffer without MgCl₂, 0.8 μL 50 mM MgCl₂ stock solution, 1 μL 10 mM PCR nucleotide mix, 0.5 μL 20 mg mL^-1 BSA (Roche Diagnostics GmbH, Mannheim, Germany), and 1 μL each of 10 μM primer and 5 ng DNA template. The thermal cycler protocol was 95°C for 5 min, 35 cycles of 95°C for 60 s, 52°C for 30 s, 72°C for 2 min and a final 7 min extension at 72°C. Amplicons were run with 1% (w/v) agarose gel to check the band size, and then sequenced at LGC Genomics GmbH (Berlin, Germany) on an Applied Biosystems^® 3730XL DNA analyzer. The 16S rRNA gene sequences of bacterial isolates were compared with the reference sequences available in the Nucleotide Database of the National Center for Biotechnology Information (NCBI) using the basic local alignment search tool (Nucleotide BLAST)¹. The accession numbers of the 68 unique strains identified at species level by the 16S rRNA gene sequences are MF664115 to MF664179 and MF677865 to MF677867 (Supplementary Table S3). Phylogenetic and molecular evolutionary analyses with the 16S rRNA gene sequences of bacterial isolates were conducted by using software MEGA 7.0 for bigger datasets (Kumar et al., 2016). The sequences were aligned by using the CLUSTALW (Thompson et al., 1997). Tree constructions were performed using the Maximum Likelihood [ML] method (Cavalli-Sforza and Edwards, 1967; Felsenstein, 1981, 1993). The robustness of the phylogenetic tree was confirmed by using 1000 bootstrap replications. Additionally, Chlorella vulgaris C-27 and Planctomycetes bacterium AS92 were included as the outgroup, and the validation of the phylogenetic neighbors was carried out by adding the 16S rRNA sequences of type strains obtained from the SILVA rRNA database project (Quast et al., 2013; Yilmaz et al., 2014). Phylogenetic tree was visualized and exported using the web-based tool Interactive Tree Of Life (iTol) (Letunic and Bork, 2011).

To test the effect of different types of carbon sources on the bacterial growth, 14 carbon sources including Alanine, Arabinose, Butyrolactam (2-Pyrrolidone), Fructose, Galactose, Glucose, Glycerol, Glycine, Lactic Acid, Putrescine, Serine, Succin Acid, Threonine, Valine were used and incorporated individually into modified OSG medium (minimal medium containing all salts required for bacterial growth) (Schnider-Keel et al., 2000) (15 μL carbon source into 135 μL modified OSG medium) at the final concentration of 0.05% (w/v). After inoculation of 5 μL bacterial culture into the 150 μL modified OSG medium containing corresponding carbon source onto a 96-well microtiter plate, followed by incubation for 2 days at 25°C, growth was measured by absorbance at 600 nm. For each microtiter plate, negative controls were included by adding LB medium without bacterial culture together with the other reagents.

The statistical analyses in this study were performed in R environment² (R Core Team, 2013). CFU values were log transformed before statistical analysis. Significant differences in log (CFU) across plant compartments, plant species, and soil types were identified using three-way analysis of variance (ANOVA) (function aov), followed by multiple comparisons (function TukeyHSD) in multcomp package (Hothorn et al., 2008). Prior to ANOVA analyses, normality of the dataset was tested by using base function sapply() modified by adding the function mystats by Kabacoff (2015), while homogeneity of variances was verified by using the function fligner() in R.

The data associated with the functional traits (biochemical assays) were subjected to a range of normalizations. First, the absorbance values obtained for bacterial fitness (abiotic and biotic resistance) and metabolic potential (growth in 14 carbon sources) were normalized by the absorbance of the bacterial growth in LB medium without other reagents at the stationary phase. Traits associated with bacterial fitness were further weighted to account for the relative abundance of each strain in the original samples. This was achieved by multiplying the absorbance data by the relative abundance of each strain, i.e., the ratio of the number of colonies identified as a specific strain to the total number of colonies obtained from the corresponding treatment. For the traits associated with the metabolic potential, the normalized data were further standardized by the maximum absorbance value observed across all carbon sources and used to calculate the niche breadth for each bacterial isolate, by summing all 14 values according to Salles et al. (2009). Regarding the traits associated with plant growth promotion, which are based on bacterial supernatant, the absorbance values were only weighted to account for the relative abundance of each strain in the original samples, as described above.

In order to show the variation of functionality along the chronosequence, the variation of individual functional traits along the chronosequence were simulated using polynominal regression performed in polynom package (Venables et al., 2016) and plotted in ggplot2 package (Wickham, 2009), and only the significant regressions were shown. The polynomial regression was used to verify the distribution of traits along the chronosequence containing five successional stages, as previously done for data referring to plant diversity and biomass, bacterial diversity and microbial activities (Schrama et al., 2013; Dini-Andreote et al., 2014, 2015; Wang et al., 2016; Salles et al., 2017). To test for correlations between bacterial growth under pH and salinity stress and the variations of soil pH and salinity along the succession, we applied Spearman’s rank-based correlational analysis. The distribution of functional traits performed by isolates from different successional stages and plant compartments was performed by principal components analysis (PCA), using the built-in R function prcomp(). Values of the individual functional traits were log transformed before PCA. The visualization of the PCA results was conducted by using the function ggbiplot() in the packages ggbiplot (Vincent, 2011).

To further verify the overall patterns of functional diversity, we generated functional diversity measure encompassing all the quantified traits, by modifications from Laliberté and Legendre (2010). A matrix containing all normalized data on individual trait for each species was used to calculate the pairwise distance between species based on Bray Curtis distance. We then calculated the average pairwise distance for each treatment (plant compartment, plant species and soil type) by pooling the treatment-specific species together. The pairwise comparisons for the overall functional diversity between two stages were tested by using Post hoc test (function posthoc.kruskal.nemenyi.test) in PMCMR package (Pohlert, 2014) after Kruskal–Wallis test.

The strength of each functional trait associated with the bacterial isolates from different plant compartments, plant species and soil types were visualized by using heatmap (function aheatmap) in NMF package (Gaujoux, 2015; Seoighe, 2015). Prior to analysis, the normalized absorbance values were further standardized by using Z-scores. Pie charts were used to summarize the proportion of bacterial isolates from one of the treatments.

The population counts on R2A agar plates of bacterial isolates from the rhizosphere soil and root endosphere of L. vulgare and A. maritima along the chronosequence were significantly influenced by soil type (F = 10.820, P = 0.002). Significant differences in population counts were respectively observed between middle (15- and 35-year stages) and late (65- and 105-year stages) successional phases for either rhizosphere or endosphere from both plants (Tukey’s HSD test, for L. vulgare, rhizosphere, P = 0.012, endosphere, P = 0.008; for A. maritima, rhizosphere, P = 0.013, endosphere, P = 0.001). In addition, significant polynomial variations were found for either compartment from both plant species (Figure 1), respectively decreasing from the initial (5-year stage) to middle phase followed by an increase toward the late phase along the chronosequence.

In order to foster variation among isolates, up to 30 colonies differing in morphology were recovered from each treatment, generating a total of 808 colonies. These were characterized by genotypic characterization using BOX-PCR, which targets repeated regions across the bacterial genome, therefore depicting differences at strain level (thus beyond 16S rRNA gene identification, Cho and Tiedje, 2000). The BOX-PCR analyses indicated that these 808 colonies could be assigned to 159 bacterial genotypes (strains) with unique BOX-PCR patterns (Supplementary Figure S1 and Table S1), which were further identified at species level by sequencing of the 16S rRNA gene, revealing a set of 68 unique bacterial species (Supplementary Table S2). From those, 59 species were unique to individual treatments and further functionally characterized using biochemical assays targeting a range of traits. The remaining 9 strains were observed in more than one plant compartment or plant species (Hafnia psychrotolerans strain CSE_16, Pseudomonas sp. CanL-3, Pseudomonas sp. 332, Pseudomonas sp. UT 6-06, Pseudomonas sp. ARCTIC-P37, Pseudomonas sp. JY-Q, Pseudomonas sp. 46 (2016), Pseudomonas fluorescens strain FW300-N2E2, and Serratia fonticola strain 51) and therefore not included to avoid biases in the analyses (Supplementary Figure S2).

According to phylogenetic analyses, the 59 unique strains were affiliated with six taxonomic groups (Supplementary Figure S3). Gammaproteobacteria was the predominant class with 45 unique strains included, from which 33.3% of the strains were isolated from the 65-year stage while the remaining ones were almost evenly distributed among other stages. Moreover, the strains belonging to Gammaproteobacteria were equally isolated from either plant species or plant compartment (Supplementary Figure S4). Actinobacteria was the second dominant phylum including 6 unique strains, from which 3 strains were obtained from 5-year stage. Among these Actinobacteria isolates, 66.7% were from A. maritima, among which 75% were from rhizosphere. The strains affiliated with Bacilli (two isolates) and Bacteroidetes (four isolates) were all isolated from the rhizosphere of A. maritima, whereas the single Alpha- and Beta-proteobacteria isolates were from the endosphere of L. vulgare. Specifically, the strains belonging to Bacteroidetes, Alpha- and Beta-proteobacteria were only isolated from 5-year stage. Thus, the diversity of bacterial isolates was highest at 5-year stage including five taxonomic groups (16 unique strains), among which Gammaproteobacteria accounted for almost half of the total amount of strains (Supplementary Figure S5). Whereas the lowest taxonomic diversity was observed at 65-year stage, where 14 out of 15 isolated strains were affiliated with Gammaproteobacteria.

When analyzed along time – as represented by the different soil stages along the primary succession – some of the individual functional traits showed significant patterns along the chronosequence. Specifically, isolates obtained from the rhizosphere showed a hump-shaped distribution for salinity stress resistance (r² = 0.323, P = 0.006) and siderophore production (r² = 0.206, P = 0.049), peaking at the 35-year stage (Figures 2B,I). The salinity stress resistance strongly correlated with soil salinity, whose variation progressively increases along the succession (Wang et al., 2016) (spearman correlation, P = 0.005; polynomial regression, r² = 0.219, P = 0.042) (Supplementary Figure S6). Similarly, hump-shaped tendencies were found for antibiotics resistance (penicillin: r² = 0.360, P = 0.003; streptomycin: r² = 0.357, P = 0.003), but those traits peaked later in the succession, at the 65-year stage (Figures 2E,F). Regarding the isolates obtained from the endosphere, we only found significant patterns for resistance to salinity (r² = 0.437, P = 0.031) and osmotic stress (r² = 0.368, P = 0.005) along the chronosequence (Figures 2B,D), showing a reverse hump-shaped pattern with a slight decrease from the initial stage to 65-year stage followed by a sharp increase toward 105-year stage, while other traits were relatively constant across the succession (Figures 2A,C,E–I). For both rhizosphere and endosphere, no significant trend of bacterial growth under different pH was observed (Supplementary Figures S7, S8).

The functional trait patterns obtained from the bacterial isolates from different plant compartments across the succession were represented using multivariate analyses (PCA, Figure 3). PC1 and PC2 accounted for 60.67% of the total variance among all traits tested in this study, largely reflecting the functionalities of bacterial isolates. PC1 was mostly represented by antibiotic resistance, and growth under different pH (especially pH = 6, 8, 9). PC2 was mainly represented by production of exoprotease and siderophore. In addition, resistance to salinity stress is associated with a PC3 (not shown), given the 45° angle with both PC1 and PC2.

The distributional pattern of functional traits of bacterial isolates from each stage was unique. Stage 105 year showed the most restricted distribution among others as most of these strains performed relatively low functionalities (Figure 3), except for two endosphere isolates showing relatively high growth level under different pH as well as resistance to salinity and osmotic stress. Conversely, the isolates from stage 35 year showed a broad functionality, given their great dispersion in the PCA. Among these small number of strains with high functional profiles, one rhizosphere isolates showed high production of exoprotease and siderophore, whereas one endosphere isolate was found to show high resistance to antibiotics and high growth level under different pH. Similarly, the dispersion level of isolates from 65-year stage was also high. However, a large number of these isolates, which are mostly endophytes, showed high growth level under different pH, while two rhizosphere isolates were found to perform high resistance to salinity and osmotic stress as well as the production of exoprotease and siderophore. The dispersion level of isolates from 5- and 15-year stages was moderate compared to others. In addition, most of the isolates these two stages showed relatively low functional profiles, except for one endosphere isolate with high antibiotic resistance and two other endosphere isolates with high resistance to salinity and osmotic stress from 15-year stage.

Among the isolates showing high biochemical properties, which were mostly from 65-year stage, Pseudomonas sp. DSM29166, Serratia plymuthica strain NBRC102599, Serratia plymuthica PRI-2C, Serratia plymuthica strain I-A-E-24, Erwinia sp. strain SH18, and Erwinia persicina B57 showed both resistance to abiotic stress and antibiotics (Figures 4, 5). Strains showing high functional profile in plant associated traits were also largely from 65-year stage (37.5%) compared to other stages, including Pseudomonas sp. R76, Pseudomonas putida AV4 and Stenotrophomonas maltophilia strain ATCC 13676 (Figure 6). Similarly, a larger proportion of the isolates showing high metabolic potential were from 65-year stage (37.5%) compared to others, including Pseudomonas sp. CH2(2014), Pseudomonas sp. 12Kp11 and Microbacterium foliorum YIM130897 (Figure 7).

To further verify the patterns of functional diversity, we generated functional diversity measures encompassing all the traits measured, allowing us to depict functional similarity among species along a continuum of soil formation as a whole (total functional diversity) or per category of function (specific functional diversity): traits associated with bacterial fitness (antibiotic and abiotic stress resistance – pH, osmotic and oxidative stress, and salinity), metabolic potential determined by growth on different carbon sources, and PGP capacity (siderophore production, exoprotease production and biofilm formation).

Analyzes of variance indicated soil type exerted significant influence on total functional diversity (P = 0.035), but not on that for bacterial isolates from either plant compartment (P > 0.05) (Table 1). In addition, the total functional diversity at the 35-year stage was higher compared with the initial and late stages, confirming the broad distribution observed in the PCA. When grouping the functional traits according to soil stages, we observed that soil type influenced the functional trait diversity associated with antibiotics resistance (P = 0.014) in rhizosphere with a significant enrichment at the 35-year stage, whereas for endophytes, none of the diversity of functional traits was affected by stage (P > 0.05) (Table 1). Plant species exerted significant effect on the functional diversity related to resistance against abiotic stress, metabolic potential and PGP capacity (respectively, P = 0.004, P = 0.001, P = 0.049) (Table 2). Moreover, the total functional diversity of isolates obtained from L. vulgare followed the development of succession, peaking at the 15-year stage (Table 1).

From the taxonomic characteristics, most of the bacterial isolates were unique to one plant species or plant compartment (Supplementary Figure S2), which coincided with the effect of plant selectivity, recruiting unique microbes in response to the specific root exudates produced by different plants (Hartmann et al., 2009; Berendsen et al., 2012) and the degree of association with plant host, given the different biochemical environment observed in the internal root tissues (Gaiero et al., 2013; Timm et al., 2015). Specifically, this selectivity was obvious for the non-dominant groups — Bacilli, Bacteroidetes, Alpha- and Betaproteobacteria, probably because of the low number of isolated species associated with these classes. We observed a predominance of Gammaproteobacteria, especially belonging to the genus Pseudomonas, in our isolates despite the utilization of general isolation media. Similar results have been found elsewhere (Marilley and Aragno, 1999; Haas and Défago, 2005; Bharathkumar et al., 2008; Hong et al., 2013), which could be explained by the important functions that the genus Pseudomonas provides for plant hosts, such as plant growth promotion (Patten and Glick, 2002; Berg and Smalla, 2009; Rashid et al., 2012) and biological control (Chin-A-Woeng et al., 2003; Garbeva et al., 2004b; Beneduzi et al., 2012). The absolute predominance of Gammaproteobacteria from 65-year stage (Supplementary Figure S5) was also confirmed in our previous study on the rhizosphere of L. vulgare where more unique OTUs belonging to Gammaproteobacteria (a total of 99) were found at 65-year stage compared to other stages (Wang et al., 2016). These results could serve as an inference of a more stabilized community structure toward the end of succession, being in accordance with the reduced community turnover in rhizosphere (Wang et al., 2016), possibly resulting from the dominance of buffering effects of soil as well as plants following the development of succession (Walker and del Moral, 2003; Dini-Andreote et al., 2014). Despite the low number of individuals belonging to other classes, it is interesting to note that the highest diversity of other taxonomic groups was observed at the initial stage, which is subjected to a constant influx of different microbes owing to marine input (Dini-Andreote et al., 2014), followed by colonization in the different plant compartments.

The isolates obtained from the rhizosphere showed significant patterns along the chronosequence, specially for functions associated with bacterial fitness, indicating their response to soil properties as well as their potential to improve plant health. Specifically, the progressively increasing pattern of salinity stress resistance, antibiotic resistance and siderophore production of bacterial isolates were in line with the soil salinity accumulation and soil nutrients enrichment following the development of succession (Dini-Andreote et al., 2015; Wang et al., 2016). Moreover, it confirms the metagenomics data generated from bulk soil samples collected from the same field, where we observed an increase in genes associated with both salinity stress and antibiotic resistance at the end of the chronosequence (Dini-Andreote et al., 2017, submitted). These traits associated with antibiotic resistance are believed to play an important role in the adaptation of microbes to complex-structured environments like the soils found at the late successional stages. Despite the concurrence between studies, it is important to notice that the patterns observed in this study are obtained for a limited number of isolates and that an increase in the survey could lead to changes in the patterns. For instance, we expect that the difference between different plant compartments and plant species would be decreased, considering that more bacteria performing similar functions would be isolated. Below we discuss the most significant results obtained by grouping the bacterial isolates according to their specific functional diversity.

Contrary to the rhizosphere isolates, the majority of functions associated with endophytes, except for resistance to salinity and osmotic stress, remained constant along the chronosequence, confirming our hypothesis that plant selectivity would be stronger in these communities. Thus, the selectivity and the buffering effect offered by plants against abiotic and biotic stresses (Hallmann et al., 1997; Rosenblueth and Martínez-Romero, 2006; Schulz and Boyle, 2006) lead to a stable community structure (Wang et al., 2016) and further functional profile for the root endophytes. Interestingly, the resistance against salinity and osmotic stress dramatically increased at the 105-year stage, which was in line with the community turnover of endophytes observed earlier when assessing L. vulgare-associated bacterial communities through molecular methods (Wang et al., 2016). This could be derived from the stress experienced by both L. vulgare and A. maritima because of the higher salinity level and severe competition with other plants at the late stage, specially Elytrigia atherica, which dominates the vegetation (Schrama et al., 2012). The reason why these two traits were synchronized in response to the soil and plant changes along the succession could be explained by the inhibitory effect of salt on plant growth, causing both ion toxicity and osmotic stress (Zhang et al., 2010; Tavakkoli et al., 2011).

By disentangling the distributional pattern of functional traits from different treatments and by grouping the functional traits into an overall functional diversity measure we observed that the functional diversity of bacterial isolates was highest at intermediate stages of succession. We could speculate that the intermediate soil nutritional level, intermediate stressful abiotic (salinity, higher flooding frequency) and biotic environment (plant competition) at the stage could lead to a selection of more versatile root-associated bacteria that can better adapt to and cope with the variety of biotic and abiotic conditions, eventually leading to an enrichment of functional diversity in the intermediate succession. This high functional diversity was, however, not associated with the taxonomic diversity, which was highest at early stages of succession. The screening of a larger number of isolates should provide a better characterization of the potential (lack of) link between taxonomy and function, which should be the subject of follow up studies. For instance, the screening of larger number of strains from the initial stage could eventually increase the functional diversity patterns, which should lead to a larger dispersion in Figure 3. However, given the low functional variation observed despite the high taxonomic diversity, we argue that is not likely to occur. Likewise, the screening a larger number of isolates should not modify the observed high functional diversity assigned to intermediate stage. The lowest functional diversity was observed in the 105-year stage, characterized by high soil nutrient content, plant biomass and low environmental fluctuations. We speculate that this low variability in abiotic and biotic stress lead to a more stable, functionally specialized microbiome. We expect that this pattern would remain the same if more isolates would be screened, given that the bacterial isolates in the dominant taxonomic group at this stage are closely related in phylogeny. Additionally, the plant species exerted an influence on the root-associated bacterial functional diversity which was higher in the A. maritima microbiome, which mainly derived from the difference of root exudates from different plants (Hartmann et al., 2009; Bednarek et al., 2010; Berendsen et al., 2012), resulting in plant specificity of microbial communities (Haichar et al., 2008; Berg and Smalla, 2009; Turner et al., 2013; Ofek-Lalzar et al., 2014).