Sampling.

To enable a robust biogeographical analysis of the patterns of bacterial colonization onTamarix aphylla leaves, we sampled leaf specimens of this tree along a 500-km transect in the southwestern United States. The collection sites were identified from previously recordedT. aphylla establishments (data provided by the participants of the Consortium of California Herbaria; ucjeps.berkeley.edu/consortium/) and selected for their relatively uniform climatic conditions (Table 1) and geographic distribution throughout the Sonoran Desert in California and Arizona (Fig. 1). The mean maximum temperatures at the time of sampling (June 2010) were in the range of 38 to 40°C.Tamarix aphylla trees are not native to the area; they were planted decades ago for decorative purposes or as windbreaks (6).

Five sites were selected and designated A to E (Fig. 1); two trees were sampled within 1 km of each other at each site. Exact GPS coordinates were recorded for each location. Branches were collected in sterile paper bags and processed within a few hours. Leaf tissue samples, 40 g each, were immersed in 40 ml sterile 0.1 M KPO₄ (pH 7), sonicated for 2 min in an ultrasonic cleaning bath, and then vortexed for 15 s. The process was repeated 6 times over a period of 30 min. The leaf wash was decanted, and bacteria were recovered on 0.45-μm cellulose nitrate filters by vacuum filtration. The filters were halved and immediately frozen for transit on dry ice.

DNA extraction.

Upon arrival at the laboratory, DNA was extracted from the frozen filters (PowerSoil kit; MoBio, Carlsbad, CA, USA). Each filter half was analyzed separately to control for contamination and PCR bias. Data from both filter halves were later combined to form one data set per sample.

Analysis of the leaf environment chemistry.

Leaf washes were also collected for chemical analysis. For each tree sample, 1 g of leaf tissue was washed and vortexed twice for 30 min in 10 ml sterile water. One leaf wash was immediately filter sterilized for electrical conductivity and pH measurements and then frozen on dry ice, while the other was decanted, filter sterilized, adjusted to pH 2 for analysis of dissolved organic carbon (DOC), and then frozen. Electrical conductivity (EC) was measured using a conductivity meter (S30 Seveneasy Conductivity; Mettler, Toledo, OH). The correlation between EC values (in mS/cm) and the concentration of Na⁺ ions in the leaf wash (in mg/liter) was established using an inductively coupled plasma optic emission spectrometer (ICP-OES) (Optima 3000; Perkin Elmer, Waltham, MA) and was found to be linear (EC = 3.4 · [Na⁺] − 10.6;R² = 0.99). pH was determined using a pH meter equipped with a combination glass electrode (model 420 ThermoOrion; Orion). DOC was determined using a FormacsHT high-temperature total organic carbon analyzer (Skalar Analytical B. V., Breda, Netherlands) following the removal of inorganic carbon by lowering the pH to <2.0.

Climatic data for June 2010 were collected from Weather Underground (wunderground.com). The data used in the analysis included average minimum temperatures and average maximum humidity. There was no recorded precipitation in any of the sites during that month.

16S rRNA gene pyrosequencing.

Details of this method have been described elsewhere (4,16,32). DNA samples were PCR amplified using primers flanking the bacterial V4-V6 hypervariable regions in the small-subunit rRNA gene. These amplicons were pyrosequenced using the Genome Sequencer FLX System (Roche) with Titanium reagents/protocols. To decrease error rates, the sequencing output was passed through several quality filters in addition to the ones built into the system (15). Reads with ambiguous nucleotides, reads lacking an exact match to the sample key or to the proximal 1064R primer sequence, and reads with more than 3 mismatches to the 565F trimming anchor were omitted from analyses.

Phylogenetic and statistical analysis.

Chimeric sequences were evaluated using UChime (USearch with the −uchime option) (7). To ensure optimal chimera identification, we ran UChime with the “reference database” option (similar to Chimera Slayer [12]) using the rRNA 16S Gold fasta file of chimera-checked reference sequences, as well as with thede novo option (default options and abskew of 5; similar to Perseus [27]). Reads identified as chimeras with either option were excluded from the subsequent operational taxonomic unit (OTU) analysis. Sequences were assigned taxonomic names using Global Alignment for Sequence Taxonomy (GAST) (16) and were subsequently separated into seven subsets for clustering according to their assignment to phylum or proteobacterial class (Alphaproteobacteria,Betaproteobacteria,Gammaprotobacteria,Firmicutes,Actinobacteria,Bacteroidetes, and all other remaining sequences). Sequences were grouped into 3% OTU clusters by the single-linkage preclustering (SLP) method (17), designed to minimize OTU inflation caused by sequencing errors, by using a pairwise alignment method and a 2% single-linkage preclustering methodology; this was followed by an average-linkage clustering based on the same pairwise alignments (17,27). The same subsets were maintained to facilitate downstream analyses.

Richness estimators were calculated using MOTHUR (31). β-Diversity was analyzed using the VEGAN package in R (26). A similarity matrix was generated using a probabilistic distance metric referred to hereinafter as the “Chao” metric (2). This metric was designed to account for unseen shared species in undersampled data that have many rare species and unequal sample sizes, typical of microbial pyrosequencing data. The validity of the measurement of β-diversity based on OTU clusters was tested by comparing this sequence-based pairwise dissimilarity to a pairwise dissimilarity calculated using the GAST taxonomic assignments. The two measurement methods were found to be nearly equivalent, with the 3% OTU-based dissimilarity approximately 10% higher than the taxonomic-based one (dissimilarity_3%OTU = 1.1 · dissimilarity_genera;R² = 0.84; see Fig. S1 in the supplemental material).

We have correlated these dissimilarity matrices with environmental parameters (Table 1) that included climatic information (temperature and humidity), leaf chemistry data (salinity, pH, and DOC), and geographic location (longitude, latitude, and altitude). Comparisons were carried out using environmentally fitted, nonmetric multidimensional scaling (NMDS) plots (R VEGAN package) (24,26). The statistical significance of the relationship between (i) species dissimilarity and geographic distance or (ii) species dissimilarity and difference in environmental parameters was tested by partial Mantel tests (R VEGAN package) (20,26). This tool, selected because the dissimilarity matrix values were not necessarily independent, was employed to isolate the effects of correlated parameters from one another.P values were calculated using 100,000 permutations on rows and columns of dissimilarity matrices.

Sequence Read Archive accession number.

All sequences have been deposited in the NCBI Sequence Read Archive (project ID SRP005695).

Bacterial diversity of theTamarix phyllosphere metacommunity.

After applying all quality filters, a total of 163,895 V4-V6 16S rRNA gene sequences, 481 nucleotides in length, have been obtained from the metacommunity of epiphytic bacteria on Sonoran DesertTamarix leaves. The reads clustered into 29,048 OTUs, based upon 3% sequence dissimilarity. Proteobacteria were the most abundant, accounting for 65,944 16S rRNA gene sequences. Of the proteobacteria, the gammaproteobacteria contributed most to the diversity, with over 10,000 OTUs assigned to this phylum. The Chao1 species richness estimator was calculated using a subset of 8,700 16S rRNA gene sequences for each site (Table 2). Estimated richness ranged between 6,862 OTUs for site B1 and 21,704 OTUs for site A1. Of the environmental parameters tested, only altitude appeared to explain some of the variation in species richness (Spearman'sS = 22.8;P = 0.004).

Taxonomic assignment of the reads identified 1,010 genera, only 22 of which were found at a relative abundance greater than 5% in at least one sample (Fig. 2; also see Fig. S3 in the supplemental material). For the sake of clarity, only these genera are included inFig. 2. The majority of these abundant genera were shared by all samples. One sample that stands out as a clear outlier is D1, which appears to be dominated by the genusBartonella, comprising over 70% of the 16S rRNA gene sequences from this site while undetected in all other locations. This organism, a known intracellular pathogen (37) of mammals, is thus very likely not a part of the native microbial phyllosphere population. As site D1 was also identified in a preliminary NMDS analysis as an extreme outlier, it was omitted from further analyses.

The Sonoran Desert phyllosphere community resembles that ofTamarix from other arid regions.

To place our results in a global context, a comparison of community composition was performed, combining the current sampling effort with previously acquired data from three differentTamarix tree species in sites by the Dead Sea, the Mediterranean, and in the Negev desert, Israel (9). As the samples from Israeli sites were obtained using a different primer set than the one used here, the community comparison was performed using taxonomic assignment (GAST) (16) and not by OTU clustering. This comparison ofTamarix phyllosphere biodiversity reveals that the bacterial community of theTamarix phyllosphere has a relatively similar composition in trees from arid regions in both Israel and the United States, as reflected by their relative branch clustering shown in inFig. 3. The U.S. Sonoran communities are clearly more similar to samples from the Dead Sea region in Israel than to the wetter “Mediterranean” trees, which include the samples from the Negev and from the U.S. island of Martha's Vineyard as outliers. Since community composition differed significantly between different climatic zones in close geographic proximity (Dead Sea versus Mediterranean), regional climate appears to play an important role in the assembly of these communities. Therefore, this global community comparison underscores the importance of our efforts to isolate the distance variable by controlling for climatic variability over our sampling transect.

Similarity in community composition depends upon both geographical distance and environmental parameters.

To visualize the similarity between sites and its correlation with environmental parameters, NMDS analysis (Fig. 4, left) was performed. The relative positioning of sites in this plot, as determined by community similarity, generally corresponded to the geographical gradient. The gradient in the different factors characterizing the sampling sites is represented here by vectors, their direction highlighting the mean direction of the gradient and their length the magnitude of correlation between the gradient and community dissimilarity. This vector fitting revealed that geographical distance indeed plays a major role, with longitude being the most significant driver of community dissimilarity. Salinity and humidity were also found to have a significant, albeit weaker, effect on community dissimilarity, even under the relatively uniform climatic conditions prevailing throughout the sampled region. Interestingly, pH, an abiotic factor known to affect bacterial diversity in soils (8,19), was the only measured abiotic factor that did not appear to be correlated with the community composition of any of the tested phyla.

Although the sampling scheme was designed to minimize environmental variability, the NMDS results nevertheless suggested some correlation between geographical distance and certain environmental variables. Indeed, the combined dissimilarity in maximum daily humidity, minimal daily temperature, and leaf salinity (calculated using the Gower dissimilarity metric) (10) was linearly correlated with geographical distance (see Fig. S2 in the supplemental material). To separate the effect of geography from the environmental factors that may be correlated with it, we subjected the data to partial Mantel tests, designed to measure the correlation between two dissimilarity matrices. The tests were performed with the purpose of establishing a partial correlation between community dissimilarity, geographical distance, and the combined differences in temperature, humidity, and salinity, using both the entire bacterial community data set and taxonomic subsets. Mantel analysis of the entire bacterial community (Fig. 5A andB) revealed no significant correlation for either geographical distance (on top of environmental conditions;R = 0.25,P = 0.17) (Fig. 5A) or environmental conditions (on top of geographical distance;R = 0.1,P = 0.32) (Fig. 5B), signifying that although the effect of geographical distance is stronger, it cannot be separated from the effect of some environmental conditions. When performing this analysis on the taxonomic subsets (adjusting for multiple comparisons), only betaproteobacterial sequences (Fig. 5C andD) exhibited a significant correlation with geographical distance, assuming uniform environmental conditions (R = 0.53,P = 0.0057) (Fig. 5C). A similar correlation was not found with environmental conditions when assuming no geographical distance (R < 0,P = 0.78) (Fig. 5D), confirming that a significant distance-decay relationship is observed for betaproteobacteria in this metacommunity. An NMDS analysis for the betaproteobacteria alone confirms this observation, with a clear east-west gradient manifested along the first axis (Fig. 4B). The trends that appear inFig. 5 can be compared to the results inFig. 1 of Martiny et al. (23), displaying the “multiple habitats and provinces” scenario for the entire bacterial population and the “one habitat, multiple provinces” scenario for betaproteobacteria only.

• 1.Bell T.2010. Experimental tests of the bacterial distance-decay relationship. ISME J.4: 1357–1365 [DOI] [PubMed] [Google Scholar]
• 2.Chao A, Chazdon RL, Colwell RK, Tsung-Jen S.2005. A new statistical approach for assessing similarity of species composition with incidence and abundance data. Ecol. Lett.8: 148–159 [Google Scholar]
• 3.Cleveland WS.1981. LOWESS: a program for smoothing scatterplots by robust locally weighted regression. Am. Stat.35: 54 [Google Scholar]
• 4.Dethlefsen L, Huse S, Sogin ML, Relman DA.2008. The pervasive effects of an antibiotic on the human gut microbiota, as revealed by deep 16S rRNA sequencing. PLoS Biol.6(11): e280doi:10.1371/journal.pbio.0060280 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 5.de Wit R, Bouvier T.2006. “Everything is everywhere, but the environment selects”; what did Baas Becking and Beijerinck really say?Environ. Microbiol.8: 755–758 [DOI] [PubMed] [Google Scholar]
• 6.Di Tomaso JM.1998. Impact, biology, and ecology of saltcedar. Tamarix spp. in the Southwestern United States. Weed Technol.12: 326–336 [Google Scholar]
• 7.Edgar RC.2010. Search and clustering orders of magnitude faster than BLAST. Bioinformatics26(19): 2460–2461doi:10.1093/bioinformatics/btq461 [DOI] [PubMed] [Google Scholar]
• 8.Fierer N, Jackson RB.2006. The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. U. S. A.103: 626–631 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 9.Finkel O, Burch AY, Lindow SE, Post AE, Belkin S.2011. Geographical location determines the population structure in phyllosphere microbial communities of a salt-excreting desert tree. Appl. Environ. Microbiol.77: 7647–7655 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 10.Gower JC.1971. A general coefficient of similarity and some of its properties. Biometrics27: 623–637 [Google Scholar]
• 11.Green J, Bohannan BJM.2006. Spatial scaling of microbial biodiversity. Trends Ecol. Evol.21: 501–507 [DOI] [PubMed] [Google Scholar]
• 12.Haas BJ, et al. 2011. Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons. Genome Res.21: 494–504 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 13.Horner-Devine MC, Lag M, Hughes JB, Bohannan BJM.2004. A taxa-area relationship for bacteria. Nature432: 750–753 [DOI] [PubMed] [Google Scholar]
• 14.Hubbell SP.2001. The unified neutral theory of biodiversity and biogeography. Princeton University Press, Princeton, NJ: [DOI] [PubMed] [Google Scholar]
• 15.Huse SM, Huber JA, Morrison HG, Sogin ML, Welch DM.2007. Accuracy and quality of massively parallel DNA pyrosequencing. Genome Biol.8(7): R143doi:10.1186/gb-2007-8-7-r143 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 16.Huse SM, et al. 2008. Exploring microbial diversity and taxonomy using SSU rRNA hypervariable tag sequencing. PLoS Genet.4: e1000255doi:10.1371/journal.pgen.1000255 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 17.Huse SM, Welch DM, Morrison HG, Sogin ML.2010. Ironing out the wrinkles in the rare biosphere through improved OTU clustering. Environ. Microbiol.12: 1889–1898 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 18.Knief C, Ramette A, Frances L, Alonso-Blanco C, Vorholt JA.2010. Site and plant species are important determinants of theMethylobacterium community composition in the plant phyllosphere. ISME J.4: 719–728 [DOI] [PubMed] [Google Scholar]
• 19.Lauber CL, Hamady M, Knight R, Fierer N.2009. Pyrosequencing-based assessment of soil pH as a predictor of soil bacterial community structure at the continental scale. Appl. Environ. Microbiol.75: 5111–5120 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 20.Legendre P, Legendre L.1998. Numerical ecology, 2nd English ed. Elsevier, Amsterdam, Netherlands [Google Scholar]
• 21.Leibold MA, et al. 2004. The metacommunity concept: a framework for multi-scale community ecology. Ecol. Lett.7: 601–613 [Google Scholar]
• 22.MacArthur RH, Wilson EO.1967. The theory of island biogeography: Princeton University Press, Princeton, NJ [Google Scholar]
• 23.Martiny JBH, et al. 2006. Microbial biogeography: putting microorganisms on the map. Nat. Rev. Microbiol.4: 102–112 [DOI] [PubMed] [Google Scholar]
• 24.Minchin PR.1987. An evaluation of relative robustness of techniques for ecological ordinations. Vegetation71: 145–156 [Google Scholar]
• 25.Necola JC, White PS.1999. The distance decay of similarity in biogeography and ecology. J. Biogeogr.26: 867–878 [Google Scholar]
• 26.Oksanen J, et al. 2008. Vegan: community ecology package. R package version 1.15-1.http://cran.r-project.org/;http://vegan.r-forge.r-project.org/ [Google Scholar]
• 27.Quince C, Lanzen A, Davenport R, Turnbaugh P.2011. Removing noise from pyrosequenced amplicons. BMC Bioinformatics12: 38doi:10.1186/1471-2105-12-38 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 28.Qvit-Raz N, Jurkevitch E, Belkin S.2008. Drop-size soda lakes: transient microbial habitats on a salt-secreting desert tree. Genetics178: 1615–1622 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 29.Qvit-Raz N, et al. 2012. Biogeographical diversity of leaf-associated microbial communities from salt-secreting Tamarix trees of the Dead Sea region. Res. Microbiol.163: 142–150 [DOI] [PubMed] [Google Scholar]
• 30.Redford AJ, Bowers RM, Knight R, Linhart Y, Fierer N.2010. The ecology of the phyllosphere: geographic and phylogenetic variability in the distribution of bacteria on tree leaves. Environ. Microbiol.12: 2885–2893 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 31.Schloss PD, et al. 2009. Introducing mothur: open source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol.75: 7537–7541 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 32.Sogin ML, et al. 2006. Microbial diversity in the deep sea and the underexplored “rare biosphere.”Proc. Natl. Acad. Sci. U. S. A.103: 12115–12120 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 33.Soininen J, McDonald N, Hillebrand H.2007. The distance decay of similarity in ecological communities. Ecography30: 3–12 [Google Scholar]
• 34.Steinitz O, Heller J, Tsoar A, Rotem D, Kadmon R.2006. Environment, dispersal and patterns of species similarity. J. Biogeogr.33: 1044–1054 [Google Scholar]
• 35.Waisel Y.1961. Ecological studies on Tamarix aphylla (L.) Karst. Plant Soil13: 356–364 [Google Scholar]
• 36.Waisel Y.1991. The glands of Tamarix aphylla: a system for salt recretion or for carbon concentration?Physiol. Plant.83: 506–510 [Google Scholar]
• 37.Zeaiter Z, Fournier PE, Ogata H, Raoult D.2002. Phylogenetic classification ofBartonella species by comparinggroEL sequences. Int. J. Syst. Evol. Microbiol.52: 165–171 [DOI] [PubMed] [Google Scholar]

Supplementary Materials