Composition of supralittoral sediments bacterial communities in a Mediterranean island

Abstract

Marine coasts represent highly dynamic ecosystems, with sandy beaches being one of the most heterogeneous. Despite the key importance of sandy beaches as transition ecosystems between sea and land, very few studies on the microbiological composition of beach sediments have been performed. To provide a first description of microbial composition of supralittoral sediments, we investigated the composition of bacterial communities of three sandy beaches, at Favignana Island, Italy, using metagenetic approaches (Terminal-Restriction Fragment Length Polymorphism, sequencing of 16S rRNA genes by Illumina-Solexa technology, functional genes detection, and quantitative Real-Time PCR). Results showed that the investigated beaches are harboring a rich bacterial diversity, mainly composed by members of classes Alphaproteobacteria, Gammaproteobacteria, Flavobacteria and Actinobacteria. The metagenetic analysis showed profiles of decreasing beta diversity and increasing richness, as well as a differentiation of communities, along the sea-to-land axis. In particular, members of Firmicutes and Proteobacteria displayed contrasting profiles of relative abundance (to decrease and to increase, respectively) along the sea-to-land axis of the beach. Finally, a search for the presence of genes related to the nitrogen and carbon biogeochemical cycle (nifH, nosZ, pmoA/amoA) detected the presence of ammonia monoxygenase sequences (amoA) only, suggesting the presence of bacterial ammonia oxidation to some extent, probably due to members of Nitrospira, but with the lack of nitrogen fixation and denitrification.

Introduction

Sandy beaches are world-wide distributed and are constantly subjected to biotic and abiotic disturbances, represented by natural and artificial bioturbation, seasonal and tidal temperature fluctuation, erosion by currents, etc. (for a review see McLachlan and Brown 2006).

In spite of their importance as an ecological transition zone between land and sea, and serious concerns for their ecological persistence (Brown and McLachlan 2002; Schlacher et al. 2007), only very recently the microbial ecology of sandy sediments, particularly concerning the submerged ones, has stirred attention (Gobet et al. 2012) and a census of web sites related to bacteria in sand has recently been published (Wackett 2013). Bacterial and fungal strains have been previously isolated from beach sediments (Teplinskaia 1978; Figueira and Barata 2007; Jin et al. 2011), but very few investigations have been done aiming to describe the taxonomic composition of these sandy beaches (Rosano-Hernandez et al. 2012), consequently hampering functional studies and biodiversity estimations in such an “extreme”, but common environment. In fact, the large fluctuations in temperature, humidity, salinity and nutrient of sandy beaches, especially in a temperate zone, could allow for consideration of these environments as nontypical, when compared to soil or water environments. Additionally, the supralittoral zone of sandy beaches may also contain human pathogens, due to human impact by recreational use of beaches, or urbanization (see for examples Bonadonna et al. 2003; Mudryk 2005; Ugolini et al. 2008).

The aim of this work is to provide a first insight into the composition of bacterial communities present in supralittoral sediments of Mediterranean sandy beaches by using three beaches at Favignana Island (Italy), which are differentially exposed to wind and water streams and to anthropic impact.

Material and methods

Sampling site description, sampling procedure and physico-chemical characteristics

Samples of subsurface sand (5 cm below the surface) were taken in summer 2011 in three beaches of Favignana Island (Italy). The Favignana Islands are part of the Protected Marine Area “Isole Egadi”, located in the Southern part of Mediterranean sea, in close proximity to Sicily Island (Italy) and is characterized by Lower Pleistocene carbonate grainstones. The three main beaches present in the Island were sampled, namely, Praja (37°55′45.62″ N, 12°19′30.66″ E), which is located in close proximity to the main urban centre of the island on the north-east coast; Lido Burrone (37°55′9.67″ N, 12°18′24.67″ E), a touristic beach on the south-west coast, and Faraglioni (37°56′42.83″ N, 12°16′46.93″ E), a beach with limited touristic use, due to lack of main roads for accessing it, on the north-west coast. From every beach, a transect from the damp band to the upper limit of the beach (thereafter referred to as the “Y-axis”) was done considering three sampling points: i) the damp band (here named as the shore-line), ii) the intermediate zone between the damp band and the upper limit of the beach (named as the mid-line) and iii) the upper limit of the beach, 1 m before the dune zone and the vegetation (named as the upper-line). For each sampling point, three samples 20-cm apart from each other, were taken along an ideal line parallel to the shore-line. A total of nine samples per beach were taken. Sampling consisted of filling completely a 50 ml polypropylene sterile tube inserted in sandy sediment from its surface down to the total length of the tube (ca. 10 cm). Due to the impossibility of immediately storing at −80°C or extracting DNA, samples were stored with an open lid at ambient temperature (25°C) for two days in the dark, trying to simulate as much as possible the natural conditions present in the beach and then stored at −80 °C prior to DNA extraction. We cannot exclude that for some samples the storage conditions could have affected the relative abundance of some taxa.

Physical-chemical characteristics are reported in Table 1. Determination of organic carbon (TOC) present in the test samples was performed following a standard protocol (Italian Official Bulletin, G.U. n ° 248 of 10.21.1999) with a Perkin Elmer elemental analyzer CHNS/O Series II model 2400. Granulometry was analyzed by particle size analysis using standard procedures (Bowles 1988; Head 1984).

Table 1 Physico-chemical features and bacterial community diversity of Favignana sandy beaches

DNA extraction, T-RFLP profiling, functional genes detection and real-time PCR

DNA was extracted from sediments, after homogenization of the total sample contained in the 50 ml polypropylene sterile tube, by using a commercial kit (Fast DNA Spin kit for soil, QBiogene, Cambridge, UK), following manufacturer’s instructions.

The 16S rRNA genes were amplified from extracted DNA of all 27 collected samples with primer pairs 27f and 1495r, as previously reported (Trabelsi et al. 2009). Terminal-Restriction Fragment Length Polymorphism (T-RFLP) procedure was then followed on purified amplification products. In particular, amplicons were digested separately with restriction enzymes TaqI and AluI and digestions were resolved by capillary electrophoresis on an ABI310 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) using LIZ 500 (Applied Biosystems) as a size standard. T-RFLP analysis was performed on two technical PCR replicates from each DNA extract, as previously reported (Mengoni et al. 2005). Only peaks present in both duplicate runs were considered for successive analyses.

Real-Time PCR for quantification of bacterial cells was performed estimating the number of 16S rRNA gene copies with bacterial primers Eub341F (5′- CCTACGGGAGGCAGCAG-3′) and Eub515R (5′-TACCGCGGCKGCTGGCA-3′) targeting the 16S rRNA gene, as reported by Simmons and coworkers (Simmons et al. 2007), using genomic DNA of Sinorhizobium meliloti Rm1021 as the concentration standard for copy number calculation (genome size 6.9 Mbp). Data obtained were compared with one-way ANOVA.

Functional diversity was evaluated by PCR amplification of nifH (encoding the nitrogenase reductase subunit), nosZ (encoding the nitrous oxide reductase gene) and pmoA/amoA (the internal fragment of particulate methane monooxygenase and ammonia monooxygenase genes, Holmes et al. 1995). For nifH the semi-nested approach described in Widmer et al. 1999 was used, following the amplification protocol reported in Giuntini et al. 2006. For nosZ, primers nosZ-F-1181 and nosZ-R-1880, were used (Rich et al. 2003) with the amplification protocol already reported (Pastorelli et al. 2011). For pmoA/amoA, primer pairs A189F/A682R and A189F/mb661R (Bourne et al. 2001; Horz et al. 2005) were used with the two-steps amplification protocol reported in Horz et al. 2005. In particular, A189F/mb661R primer pair specifically amplify pmoA sequences and not amoA sequences, allowing for discrimination between methane monooxygenase and ammonia monooxygenase genes (Bourne et al. 2001). Genomic DNA of Sinorhizobium meliloti was used as a positive control for nifH and nosZ amplification, while for pmoA/amoA genomic DNA of Nitrosomonas europaea was included as a positive control in PCR amplification.

Metagenetic analysis of 16S rRNA gene amplicons

DNA aliquots extracted from Faraglioni beach samples were pooled together, with respect to the position along the sea-to-land transect to obtain three samples, each composed by the triplicate samples of shore-line, mid-line and upper-line samples, respectively. The variable V3 region of the 16S rRNA gene pool of total bacterial community was amplified from each sample DNA with primer pairs V3-338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and V3-533R (5′- TTACCGCGGCTGCTGGCAC-3′) (Huse et al. 2008). PCR conditions were as described in Sogin and coworkers (Sogin et al. 2006). Ten independent PCR reactions per sample were done, then pooled together to produce three representative PCR amplicon libraries for each environmental DNA sampling point (shore-line, mid-line, upper-line). Products were resolved by agarose gel electrophoresis and bands were purified with a MinElute Gel Extraction Kit (Qiagen, Inc.). Quality and quantity of products was assessed on a Biophotometer (Eppendorf).

Massive sequencing was performed by Illumina-Solexa technology (Gloor et al. 2010; Bartram et al. 2011) with the pair-end protocol on an Illumina HiSeq2000 machine by Beijing Genome Institute sequencing service (www.genomics.cn/). Sequences are deposited in the Bioproject database (http://www.ncbi.nlm.nih.gov/bioproject/) with ID: 234346.

Statistical analyses and processing of T-RFLP data

Analysis of T-RFLP profiles was performed as previously reported (Pini et al. 2012). Statistical analyses were performed on a binary matrix obtained by linearly combining data from the two restriction enzymes as previously reported (Mengoni et al. 2009). Ribotypic diversity, as the number of Terminal Restriction Fragments (T-RFs) identified in each sample, was used as an estimator of richness as reported previously (Mengoni et al. 2009). Cluster (UPGMA) analysis and Canonical Correlation Analysis (CCA) were done with Past software (Hammer et al. 2001) on the Jaccard similarity matrix obtained from binary T-RFLP profiles. To test the distribution of the variance of T-RFLP profiles within and among beaches and on the Y-axis, AMOVA (Analysis of Molecular Variance, Excoffier et al. 1992) was applied using Arlequin 3.11 software (Excoffier et al. 2007). AMOVA was used as a statistical methodology alternative to the classical analysis of variance (ANOVA). AMOVA is more flexible for biological data than classical ANOVA because it does not require a prior assumption of normality of the dataset and statistical significance is computed through a permutation test (Excoffier et al. 1992; Mengoni and Bazzicalupo 2002), and it has been widely applied to analyze community profiles in molecular microbial ecology (for examples see Dalmastri et al. 1999; Tabacchioni et al. 2000; Mengoni and Bazzicalupo 2002; Mengoni et al. 2009; Pini et al. 2012).

Sequence analysis and bioinformatics

Sequence reads (of approx. 100 bp in length, covering the V3 region of 16S rRNA gene) were subjected to a first quality control step in order to eliminate low complexity reads and low quality bases. Three other quality control steps were performed: i) firstly, the presence of Illumina adaptors and primers was checked by using a standard procedure, as reported in the FASTQC manual (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) that is by collecting sequences of 50 bp that were overrepresented in the reads file and analyzing them using blast against the nt database; ii) a second step was performed using a dynamic trimming algorithm to trim the poor quality bases in all the reads of the samples, in order to delete possible ambiguous bases. A quality cut off of 28 was used to obtain an error percentage at least of 0.16 %; and iii) finally, sequences files were analyzed using FastQC in order to check for increased quality of reads.

For sequence analysis and taxonomy classification, reads with less than 50 bp were deleted from the datasets. Then the resulting reads, with length > 50 bp were analyzed by using the RDP Classifier and taxonomically assigned. The results obtained by the RDP Classifier were imported into MEGAN (MEta Genome Analyzer, Huson et al. 2007) to calculate a taxonomic classification of the reads. The taxonomic classification obtained was collapsed to different taxonomic levels (phylum, class, order, family and genus) with the purpose of analyze the absolute read abundance attributed to that taxonomic level on each dataset. Recommendations for thresholds as in Mizrahi-Man et al. 2013 were followed, which allow an error rate up to 5% at the genus level using a confidence threshold of 95% (Mizrahi-Man et al. 2013).

To assess taxa richness, rarefaction analyses were performed using the R package Vegan [http://cc.oulu.fi/~jarioksa/softhelp/vegan.html], after collapsing reads to the taxonomic levels of genus, family, order and class, with a probability of assignment of 80% or above.

Finally, data obtained from MEGAN were analyzed using R (package vegan - http://cran.r-project.org/web/packages/vegan/index.html). Firstly, absolute abundances were transformed to relative ones using this formula:

$$ {X}_{ij}^{norm}=\raisebox{1ex}{${X}_{ij}$}\!\left/ \!\raisebox{-1ex}{${X}_{j+}$}\right. $$

Where: i and j are the matrix rows and columns, respectively, X _ij is the value in the row i and the column j, X ^norm _ij is the normalized value in the row i and the column j and X _j+ is the sum of all value in the column j. Then, a heat map was made for every different taxonomic level by using custom R scripts (available upon request) and the graphic package ggplot2 (http://cran.r-project.org/web/packages/ggplot2/index.html). Relative abundances were then used to make a clustering using the UPGMA algorithm and “Bray-Curtis” distance. At phylum level the absolute abundance was resampled in order to check if the differences between the datasets obtained were statistically significant. A number of 10,000 resamples was done for every pair of datasets. Then, the differences obtained with the 10,000 resamples were tested against the real difference between each samples pair using a Wilcoxon signed rank test with a 95% confidence interval (a p-value less than 0.05 was considered significant) (Figure S1). For each phylum a boxplot of the 10,000 differences obtained was plotted for each samples pair (3 boxplots in total) using R graphics package ggplot2 (see above) (Figure S2).

As previously suggested (Hill et al. 2003; Schloss and Handelsman 2006; Shaw et al. 2008; Gerber et al. 2012) there is not a unique indicator of community diversity, consequently four diversity indices were then calculated (Hill et al. 2003): inverse Simpson index (defined as D⁻¹), Shannon-Weaver index, Richness index (as estimators of alpha diversity), Evenness index and beta-diversity (defined as “total number of taxa in one site”/”means of taxa count in all sites” -1).

All the script and classes used in this work were written in Java or R and are available upon request.

Results

T-RFLP profiling of bacterial community function and diversity

Application of 16S rRNA gene T-RFLP bacterial community profiling to the 27 total DNA samples (nine samples for three sampling points per locality) allowed the identification of 30 different T-RFs (Terminal-Restriction Fragments), two of which were present in all samples, while the others were detected in 1 –11 samples. The ribotypic diversity of communities (as number of T-RFs, Table 1) was relatively low and varied from 4.7 mean T-RFs (Praja upper-line) to 9.7 (Faraglioni shore-line). In general shore-line samples had more T-RFs than upper-line samples. Bacteria titres (no. of cells/g of sand), estimated by quantitative PCR, were not different (one-way ANOVA) between samples and varied from 1.3 ± 0.6∙10⁵ cells/g (Faraglioni shore-line) to 1.2 ± 1.1∙10⁴ cells/g (Lido Burrone upper-line). No relationships between organic carbon content, ribotype richness and bacterial cell estimates (as 16S rRNA gene copies) were found (data not shown).

CCA (Fig. 1a) showed that samples are mainly grouped according to their position along the Y-axis (from left to right) and that humidity is clearly related with such a axis, while organic carbon content is more related to differences among localities. However, also some groupings related to the different localities (beaches) were detected, in particular when cluster analysis was applied (Fig. 1b). In fact, UPGMA dendrogram highlighted a strong differentiation of Faraglioni beach samples from Lido Burrone and Praja samples. In particular, one main cluster including samples from Praja and Lido Burrone only, while another one including all samples from Faraglioni and some samples from Lido Burrone, were identified. In agreement with CCA, in this latter cluster, all upper-line samples from Faraglioni and Lido Burrone grouped together at the Jaccard similarity >0.50. Finally, it is interesting to notice, that in both CCA and UPGMA, several triplicate samples from the same point along the Y-axis clustered together at theJaccard similarity > 0.60 (i.e., Faraglioni mid-line, Lido Burrone mid-line, Praja upper-line, Praja shore-line, and Faraglioni sore-line). To quantitatively evaluate the contribution of Y-axis position and of locality (single beach) an AMOVA was carried out on T-RFLP profiles. A two-levels hierarchical partition (among sampling points and within sampling points partitions) showed that most of the variance is due to differences between sampling points (73.37%, P < 0.0001). On a three-levels hierarchical analysis (Table 2), localities and Y-axis similarly contributed to bacterial community differentiation (8.9 and 8.5% of variance, respectively).

Fig. 1

Pattern of ordination of 16S rRNA bacterial community T-RFLP profiles. a Nonmetric Muldimensional Scaling (nMDS); b) UPGMA clustering. Cophenetic correlation coefficient = 0.9098. In CCA each point represents a DNA sample for the beaches of Faraglioni, Praja and Lido Burrone; squares, mid-line samples; circles, shore-line samples; triangles, upper-line samples. Black-filled symbols, Faraglioni, grey-filled symbols, Lido Burrone; white-filled symbols, Praja. Vectors for humidity (%) and organic carbon (%) are reported

Table 2 Analysis of Molecular Variance (AMOVA) of T-RFLP profiles for three localities (Lido Burrone, Praja, Faraglioni) and three environments for locality (shore-, mid-, upper-line)

To evaluate the functionality of bacterial communities in terms of performance, some of the steps of the biogeochemical cycle of nitrogen and of carbon were checked. In particular, three gene fragments, for which “universal” primers have been developed, were chosen: i) nifH (encoding the nitrogenase reductase subunit), ii) nosZ (encoding the nitrous oxide reductase gene), and iii) pmoA/amoA (the conserved internal fragment of particulate methane monooxygenase and ammonia monooxygenase genes). Amplicons were detected only for pmoA/amoA with primer pair A189F/A682R in Faraglioni and Praja shore line samples. Semi-nested reamplification with A189F/mb661 primer pair did not produced amplification products, suggesting that most of first amplification amplicons were due to amoA-related sequences.

Diversity of 16S rRNA gene amplicon libraries from Faraglioni beach

Since T-RFLP analysis showed a pattern of community diversity along the Y-axis, we focused on one locality only for an in-deep taxonomic investigation by massive sequencing of 16S rRNA gene amplicons. We selected the locality of Faraglioni because of the lower anthropic impact and higher organic carbon present, with respect to the other two localities. A total of 380,080, 391,008, and 382,584 reads for shore-line, mid-line and upper-line samples, respectively, were analyzed. Results (Table 3, Fig. S3) indicated differences in the coverage of diversity of samples along the Y-axis. In particular, upper-line was the most diverse, while mid-line and shore-line samples were less diverse. In order to have an idea of the number of unseen taxa present in the samples, the nonparametric estimator Chao I was also calculated (Table 3), allowing to show an underestimation for all the three samples along the Y-axis, which accounted for ∼44 Families. Considering higher taxonomic levels (order and class), the difference between observed and predicted values of richness were clearly reduced. In general, Richness and inverse Simpson indices were in agreement with the above mentioned coverage data, showing a profile of increasing diversity along the Y-axis from shore-line to upper-line, while Evenness and Shannon-Weaver indices have a similar profile and did not allow for resolving the variability on the Y-axis.

Table 3 Diversity indices of samples from Faraglioni beach after collapsing 16S rRNA gene sequences to different taxonomic levels

Bacterial phyla composition and sea-to-land differences of Faraglioni beach

The pattern of relative phyla abundance along the Y-axis are reported as a heat map in Fig. 3. The clustering of the heat map patterns indicated that at the phylum level upper-line and mid-line are more similar to each other, than to shore-line. The most abundant phyla were Proteobacteria, Actinobacteria and Bacteroidetes (Fig. 2). The most abundant classes were Alphaproteobacteria and Gammaproteobacteria, then followed by Flavobacteria and Actinobacteria (Fig. 2). Then, at the order level, Flavobacteriales, Actynomicetales, Rhizobiales, Rhodobacterales, Bacillales, Alteromonadales, Chromatiales, and Oceanospirillales were the most represented (Fig. 2).

Fig. 2

Cluster analysis of the taxonomic composition of bacterial communities of the three samples of Faraglioni beach. Clustering was performed on bacterial relative abundance patterns by applying the UPGMA algorithm. The patterns of relative abundance were used to generate a heat-map and the occurrence values are reported as a barchart upon the heat-map. a Phylum level; b) Class level; c) Order level

Different trends of abundance of phyla were detected along the Y-axis (Fig. 3). In particular, Firmicutes and Proteobacteria presented increasing and decreasing trends, respectively (Wilcoxon test p-value < 0.05), while Actinobacteria showed lower abundance in the mid-line compared to the other two samples (Wilcoxon test p-value < 0.05). Bacteroidetes did not show differences. The two phyla of Fusobacteria and Nitrospirae were excluded from the analysis, since they have too few hits.

Fig. 3

Identification of abundance trends in bacterial Phyla along the Y-axis. Panels show the trends of abundance of main Phyla along the Y-axis. Shown Phyla have observed differences that are significant at P < 0.05 after a resampling analysis and Wilcoxon test. For each plot abscissa is indicated the Y-axis sample, while on the ordinate is the number of reads accounting for that phylum

Variability of Proteobacteria and Firmicutes along the Y-axis of Faraglioni beach

Since Firmicutes and Proteobacteria showed a contrasting profile, a more detailed taxonomic analysis was performed in these phyla down to the family level. For Firmicutes, most of the hits were due to members of the spore-forming aerobic bacteria of the class Bacilli, order Bacillales (Fig. 2) and family Bacillaceae (Fig. 4), which were abundant in the shore-line sample. For Proteobacteria, shore-line and mid-line samples were more similar toeach other, contrarily to the pattern observed for the overall phyla composition and for Firmicutes. Members of family Rhodobacteraceae were the most abundant and did not show high variability along the Y-axis. Other less abundant families were Alteromonadaceae, Rhodobacteraceae, Erythrobacteraceae, and Ectothiorhodospiraceae (Fig. 4).

Fig. 4

Cluster analysis for Firmicutes and Proteobacteria. The analysis was performed as those reported in the legend of Figure 3, by inspecting the diversity at family level for Firmicutes (a) and Proteobacteria (b)

Finally, to allow one to infer additionally potential functional activities related to the bacterial families detected, the most abundant families for Firmicutes and Proteobacteria, namely Bacillaceae, Alteromonadaceae, Rhodobacteraceae, Erythrobacteraceae, and Ectothiorhodospiraceae were investigated collapsing reads at genus level. Observed relative genera abundances are reported in Fig. 5. For Bacillaceae, two main genera were detected, Bacillus and Halobacillus. For Alteromonadaceae, which were mainly found in shore-line samples and were practically absent in mid-line and upper-line samples, genus Marinobacter was the most abundant. Within Rhodobacteraceae, ten main genera were represented, with Sulfitobacter being the most abundant in upper-line samples. For Erythrobacteraceae, the genus Erythrobacter was the most represented, while for Ectothiorhodospiraceae, Thiohalospira was the most abundant genus.

Fig. 5

Composition of genera of the main families of Alphaproteobacteria and Bacillales. The relative proportion of genera for Bacillaceae, Alteromonadaceae, Rhodobacteraceae, Erythrobacteraceae, and Ectothiorhodospiraceae is reported

Discussion

The supralittoral belt is a highly dynamic transition zone between sea and land and is characterized by a sharp passage from a humid zone (the damp band) to an arid area, subjected to strong fluctuations of salinity and temperature. This quite narrow (often a few meters in the Mediterranean area) zone is relatively neglected by biologists as a reservoir of biodiversity (in comparison to the dune ecosystem). On the contrary, beach sediments (especially those devoted to recreational use) are usually subjected to government controls, concerning the environmental quality and the presence of pathogens. Despite the fact that in sensitive areas, as Marine Protected Areas, all the interventions on littoral, sublittoral, dunal and retrodunal environments are usually performed after a careful check of protected species, concerning the supralittoral zone, most parts of nourishment are performed with no or little attention for micro-, meio- and endofauna. However, interventions (as beach cleaning, sand movements with caterpillars, etc.) are normal practices, especially in touristic areas (McLachlan and Brown 2006; Schlacher et al. 2007). Nevertheless, sandy beaches are an ecosystem linked to the degradation of organic material of terrestrial and marine origin (Misic and Fabiano 2005; Schlacher et al. 2008; Ugolini and Ungherese 2012). Surprisingly, here we have shown that, despite the low organic carbon present (for reference with agricultural soils, see, for instance, Ludwig et al. 2011) and the extreme conditions, the sampled sandy beaches contain appreciable number of bacterial cells (from 10⁴ to 10⁵ cells/g as detected by the qPCR approach), lower than soil but comparable to those found in the plant endosphere (Mengoni et al. 2005, 2009; Pini et al. 2012) and lower than those found in La Jolla (California) sandy beaches (Glavin et al. 2004) and in Mediterranean sea water (Zweifel and Hagstrom 1995). Moreover, a high bacterial diversity was also found, suggesting that bacterial communities could have active roles in biogeochemical cycling on a sandy beach ecosystem. The detected biodiversity is mainly ascribed to Proteobacteria, Actinobacteria and Bacteroidetes phyla and to the orders of Flavobacteriales, Actynomicetales, Rhizobiales, Rhodobacterales, Bacillales, Alteromonadales, Chromatiales, and Oceanospirillales. Within these groups, several well known marine taxa (Rusch et al. 2007) can be found (e.g., in Oceanospirillales, Chromatiales and Rhodobacterales), suggesting that a consistent part of sandy beaches microbiota is of marine origin.

Moreover, by comparing bacterial community fingerprinting (T-RFLP profiles) of sandy beaches of three different localities we have shown a considerable heterogeneity of samples, either due to location and to the position along the Y-axis. In particular, locations showed different organic carbon content (see Table 1), which could be related to the differences observed between communities. However, no relationships between organic carbon percentage and ribotype richness was observed.

For what is of concern for the differentiation of communities along the Y-axis, the metagenetic analysis of Faraglioni beach showed a decreasing beta diversity and an increasing richness (alpha diversity) along the sea-to-land axis, suggesting a tendency i) for a lower differentiation when proceeding along the axis (may be due to more homogeneous environmental parameters), and ii) for an increasing number of taxa (may be linked to proximity with soil). Among the identified bacterial groups with differences of abundance along the axis, two main bacterial phyla have opposite trends of abundance along the Y-axis: Firmicutes and Proteobacteria. Within Firmicutes, the main portion was constituted by sequences attributed to the spore-forming aerobic taxon of Bacillaceae, more abundant in the damp band, with genera Bacillus and Halobacillus being the most represented. Concerning Proteobacteria, the alphaproteobacterial family Rhodobacteraceae was the most abundant along all the transect. For this family a large number of sequences from the marine sulfur-oxidizing genus Sulfitobacter (Pukall et al. 1999) was detected, allowing one to hypothesize that some sulfur cycle may also occur on supralittoral sediments. Actually, most of bacterial taxa known to be abundant in sea water were also well present in the supralittoral sediments, as Erythrobacteraceae, Ectothiorhodospiraceae and Alteromonadaceae. Erythrobacteraceae (as also Rhodobacteraceae) is a family of Alphaproteobacteria composed by strains mainly isolated from aquatic environments, and the genus most represented in our dataset was Erythrobacter (see for instances (Lee et al. 2005 and references therein). Ectothiorhodospiraceae is a group of halophilic and haloalkaliphilic Gammaproteobacteria (Tourova et al. 2007), which includes strains with known “extremophyilic” phenotypes but also able to fix nitrogen and participate in sulfur and iron biogeochemical cycles (see for instance Hallberg et al. 2011). Indeed, the most abundant genus belonging to Ectothiorhodospiraceae found in our dataset was Thiohalospira, a genus of recently discovered chemolithoautotrophic, halophilic, sulfur-oxidizing gammaproteobacteria (Sorokin et al. 2008). Alteromonadaceae is another well known gammaproteobacterial family of marine strains (Ivanova and Mikhailov 2001). Here the genus Marinobacter, which include strains able to degrade hydrocarbons (Cui et al. 2008), was indeed dominant. Interestingly, within the most abundant phyla, also Acidobacteria and Actinobacteria were represented, with a discontinuous pattern along the Y-axis. In particular, Actinobacteria were less abundant in the mid-line sample, possibly implying the presence of both marine taxa for this phylum (Bull et al. 2005) (mainly colonizing the shore-line sample) and soil taxa (mainly colonizing the upper-line metasample). Intriguingly, Bacteroidetes, were abundant and equally well represented along the Y-axis samples, implying that bacteria in this phylum are not influences by the strongly different environmental conditions encountered along the axis (e.g., water availability etc.). It is worth noting that anaerobic Bacteroidetes are a main constituent of human gut flora, but recently members of this group have been found in marine sediments (Green-Garcìa and Engel 2012). Among the other bacterial phyla, Cyanobacteria and Planctomycetes were mainly present in the mid-line sample, suggesting the mid-line as a potential challenging environment in which primary colonization (for instance by cyanobacteria) may occur (Gorbushina and Broughton 2009) and on which the peculiar cellular differentiation of Planctomycetes in a sessile and in a pelagic form, may help colonization (Fuerst, and Sagulenko 2011). Actually, Planctomycetes have been found as dominant in some intertidal sediment communities (Musat et al. 2006). Other phyla, as Chloroflexi, Deinococcus − Thermus, Spirochaetes, Verrucomicrobia and Clamydiae also show differential presence along the beach Y-axis, and may represent peculiar adaptation or marker taxa for the upper-line (e.g., Clamydiae and Verrucomicrobia). These phyla are relatively rare in our metagenetic dataset but, as recently shown for the rhizosphere (Dohrmann et al. 2013) as well as for marine sediments (Gobet et al. 2012), rare taxa may better explain ecological adaptation of communities than abundant taxa and may provide a “seed bank” for bacterial colonization.

Finally, concerning the possible functions promoted by bacterial communities here investigated, we did not detect the presence of genes related to nitrogen fixation (nifH), denitrification (nosZ) and methane oxidation (pmoA). However, some samples showed the presence of ammonia monoxygenase (amoA) genes. By comparison with the metagenetic sequencing of Faraglioni beach, we cannot a priori exclude that the (very few) detected Nitrospira may be mainly responsible for amoA genes. Nitrospira is a well known taxon of marine ammonia-oxidizing bacteria (see for instance Haaijer et al. 2013 and references therein). We can then speculate that bacterial ammonia oxidation could be occasionally present in sandy beaches and may contribute to the input of nitrite in the sandy beach ecosystem. However, more sampling along the seasons, as well as dedicated molecular studies with metagenetic sequencing of amoA genes, will be needed to infer possible functionality in the nitrogen biogeochemical cycle of sandy beach bacterial communities.