Taxonomic and functional analyses reveal existence of virulence and antibiotic resistance genes in beach sand bacterial populations.

Coastal sands are important natural recreational facilities that have become hotspots for tourism and economic development. However, these sands harbour diverse microbial assemblages that play a critical role in the balance between public health and ecology. In this study, targeted high-throughput sequencing analysis was used to identify sand-borne bacterial populations at four public beaches in Durban. The effect of heavy metal in shaping the distribution of bacterial metacommunities was determined using canonical correspondence analysis (CCA), while the functional gene profiles were predicted using PICRUSt2 analysis. Sequences matching those of the bacterial phylum Proteobacteria were the most abundant in all samples, followed by those of the phyla Firmicutes, Actinobacteria, Bacteroidetes, and Gemmatimonadetes. Genus-level taxonomic analysis showed the presence of 1163 bacterial genera in all samples combined. The distribution of bacterial communities was shaped by heavy metal concentrations, with the distribution of Flavobacteria, Bacteroidia, and Deltaproteobacteria influenced by Pb and Zn, while B and Cr influenced the distribution of Clostridia and Gammaproteobacteria, respectively. Identified antibiotic resistance genes included the peptidoglycan biosynthesis gene II, III, IV, and V, as well as the polymyxin resistance gene, while the virulence genes included the sitA, fimB, aerobactin synthase, and pilL gene. Our findings demonstrate that beach sand-borne bacteria are reservoirs of virulence and antibiotic resistance genes. Contamination of beach sands with heavy metals selects for both heavy metal resistance and antibiotic resistance in beach sand bacterial communities. Children and immunocompromised people engaging in recreational activities on beaches may be exposed to higher risk of infection.

Introduction

Sandy beaches are of prime importance for human recreation, tourism, and the development of coastal economic zones (Jonah et al. 2015). Worldwide, coastlines with long stretches of clean, and sandy beaches have become major economic zones, with tourist expenditures on accommodation, food and drink, entertainment, and other services and goods topping US$1260 billion each year (UNWTO 2016). Pre-COVID-19 predictions had placed global coastal tourism associated expenditures at above US$5 billion a day by the year 2020 (Orams 2003), underpinning the importance of sandy beaches to world economies. More people use sandy beaches than any other type of seashore as they provide the most productive fishing grounds, and offer perfect opportunities for sand bathing, a practice most common among the young and old beachgoers alike. However, sandy beaches are not just piles of sand; they harbour their own micro-ecosystems. They receive large inputs of organic matter supplied by the seawater, consisting of phytobenthos assimilates, and products washed and leached out from seaweeds, animal faeces, and remains of plants and animals (Mudryk et al. 2013). This creates optimal conditions for the growth of a high population of organisms such as small invertebrates, bacteria, actinomycetes, fungi, yeast, virus, algae, and diatoms (Zakaria et al. 2011; Whitman et al. 2014; Di Piazza et al. 2017), making beach sand a potential reservoir for aetiological agents of disease (Solo-Gabriele et al. 2015). Since sandy beaches are dynamic and sensitive places where life is under pressure, these sandy beach dwelling microorganisms exhibit remarkable physiological and behavioural adaptation to changing environmental conditions (Defeo et al. 2009). Microbiological assemblages in sandy shorelines play crucial roles in the decomposition of organic matter and pollutants, and nutrient mineralization and recycling, ecosystem services which render sandy beaches safe for beachgoers (Amaral et al. 2016). At the same time, coastal sediments are the ultimate sink of heavy metals discharged into aquatic environments, and their analysis is thought to offer a more convenient and more accurate means of determining the degree of pollution (Al-edresy et al. 2019). This is because heavy metal toxicity has a double effect on environmental microbial assemblages. First, heavy metal toxicity selects for heavy metal tolerance/resistance (Dickinson et al. 2019). Second, it selects for antibiotic resistance in microbial pathogens (Sabry et al. 1997) by two main mechanisms, namely co-resistance and cross-resistance (Nguyen et al. 2019). The microbial risk posed to sand bathers is therefore heightened if such microbes harbour antibiotic resistance genes.

The establishment of major urban centers in close proximity to the majority of sandy beaches has translated into increased anthropogenic pressure upon these natural resources, either as a result of discharge of sewage effluents and storm runoff, or as a direct consequence of recreation associated pollutions such as bather faecal shedding (Orams 2003). Also, due to tidal wave action, the microbiological quality of beach sand is positively correlated to the microbiological quality of the beach water (Weiskerger et al. 2019). Where a beach is located close to a sewage treatment plant outfall, contact with beach sand increases the risk of gastrointestinal illness by a factor almost similar to that of coming into direct contact with sewage polluted water (Devine 2014). Beachgoers are impacted by contaminated beach sands either indirectly by degrading beach water quality through cycles of deposition and resuspension of pathogens between sand and water, or more directly through physical contact with/or ingestion of sands (Halliday et al. 2014). Sensitive populations such as children, the elderly, or those with a weakened immune system are particularly at risk for long-term effects. From a public health perspective, knowledge of the microbial assemblages inhabiting recreational sand beaches could lead to the evaluation of the levels and trends of contaminants, as well as following human contact with sand; to an assessment of the effects on public health. Currently, beach advisories and closures are issued depending on faecal indicator bacterial densities in the beach water column (Zhang et al. 2015), with relatively less attention paid to their densities in beach sands. However, at times bacterial counts in beach waters can fall to levels considered safe for swimming while higher densities remain in adjacent beach sands, exposing sand bathers to heightened risk compared to swimmers (Whitman and Nevers 2003). Therefore, this study was aimed at assessing the bacterial communities in the sands of four popular recreational beaches which form the interface between the shoreline of the city of Durban and the Indian Ocean. A 16S rDNA-targeted high-throughput sequencing approach was used to determine the bacterial community structure and composition across the four beaches. Furthermore, the study determined the possible correlations between heavy metal concentrations at the various study sites and the distribution of bacterial metacommunities. Finally, the study predicted the possible functional gene profiles of the bacterial metacommunities using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) analysis. Though such microbial communities are highly dynamic in nature, such that a study carried a year ago may not necessarily be used to predict the microbial quality of recreational beach sand today, no study has ever been carried out to determine the bacterial and fungal communities of Durban’s sand beaches. This is a significant knowledge gap given the popularity of Durban’s beaches; a city that attracts the 6th most number of tourists into the African continent due to its year-round friendly climate and long expanses of sandy beaches. To the best of our knowledge, this study stands to be the first ever to characterise the bacterial communities inhabiting the sandy beaches of Durban using high-throughput sequencing analysis.

Methodology

Description of study site and sample collection

Seventy-nine percent of the South African coastline is composed of sandy beaches. The city of Durban boasts a wealth of coastal resources along its 97 km of coastline with the Indian Ocean. These include rocky shores, mangrove forests, coral reefs, coastal forests, wetlands, and sandy beaches. Of these, the least studied and most underappreciated are sandy beaches. Durban beaches are tourism hotspots owing to good weather that the City experiences all year-round, with average winter temperatures ranging from 11 °C to 24 °C, while summer temperatures average between 20 °C and 29 °C. Durban’s beaches are characterised by vast expanses of white sands, and gentle slopping coastlines, making them convenient for sand bathing, surfing, and swimming. Tourism is very important to the local economy. Samples were collected from the South Beach, Harbour Beach, Central Beach, and North Beach (Fig. 1). The later three beaches are in the range of approximately 3 km from each other, while South Beach is furthest from the rest, approximately 15 km. Approximately 100 g of sand sample were collected at each beach from five points about 5 m apart. A core sampler was used to collect sand samples in a vertical fashion from the surface down to about 10 cm deep at the foreshore zone. Sand samples from each beach were put in specimen bottles, labelled according to the collection sites, and transported to the laboratory at the University of South Africa in cooler boxes containing ice at 4 ºC for further processing.

Fig. 1

Map showing the location of the four beaches from which samples were collected, along the Durban Indian Ocean coastline, South Africa

Analysis of sand heavy metal concentration

Sand heavy metal concentrations were analysed following a previously described method (Sibanda et al. 2019). Briefly, the sand samples were dried at 105 °C for 2 h and then placed in a desiccator for 18 h after which ≈ 0.5 g of dried samples were each weighed and transferred into separate Mars6 microwave digestion vessels. To each, 9 mL analytic grade concentrated HNO3 and 3 mL concentrated HCl were added, after which the samples were heated to 175 °C for 20 min, and holding at that temperature for another 10 min. The samples cooled, filtered, and transferred into 50 mL volumetric flasks. Sample volumes were topped up with double-distilled water. Following this, 20 mL of the liquid samples were filtered and further acidified by adding 200 uL of concentrated HNO3. A multi-elemental SRM was used to prepare calibration standards, which ranged from 0.1 to 10 ug/L, and included a calibration blank. Analysis was performed on a PerkinElmer NexION 300X Q-ICP-MS. After start-up and warm-up, the instrument was auto-tuned to maximise sensitivity and minimise double charge and oxide interference. Calibration standards and samples were measured under KED (Kinetic Energy Discrimination) mode. The software generates a calibration curve after measuring solutions of known concentrations (standards), against which the unknown samples are measured.

DNA extraction and polymerase chain reaction amplification

For DNA extraction, about 10 g of each beach sand sample was initially mixed with 10 ml phosphate-buffered saline (PBS, pH 7.4). The mixtures were agitated by vortexing and allowed to stand for an hour at room temperature to dislodge bacterial cells adhering to sand particles. Following this, 400 µl supernatant aliquots were then used as samples in the extraction of total genomic DNA using Faecal/Soil Total DNA™ extraction kit (Zymo Research Corporation, CA, USA) according to the manufacturer’s instructions. The extracted DNA was first amplified using the universal bacterial 16S rDNA primers (27F and 1492R) to cover the whole variable region under the following PCR conditions: initial denaturation at 95 °C for 5 min, followed by 32 cycles of melting at 95 °C for 1 min, annealing at 55 °C for 1 min, and elongation at 72 °C for 1 min. The last step was final product elongation at 72 °C for 10 min, followed by cooling at 4 °C. Subsequently, a second PCR run was carried out using the 27F and 518R primer sets, with overhanging adapter sequences that are compatible with Illumina index as described by Selvarajan et al. (2018).

Library preparation and sequencing

The resultant PCR products were cleaned and concentrated using AMPure XP beads (Beckman Coulter, Agencourt Bioscience Corporation, Massachusetts, USA) according to the manufacturer’s instructions. Following the purification step, Illumina sequencing adapters and dual-index barcodes were added to the amplicon targets using a full complement of Nextera XT indices (Illumina, Inc. San Diego, CA, USA) through 8 cycle PCR as follows: 95 °C for 3 min, 8 cycles of {95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s}, final extension at 72 °C for 5 min, and cooling at 4 °C. The resulting PCR product was cleaned again as already described. Fragments with an approximate size of 630 bp were validated using a Bioanalyzer DNA 1000 chip (Agilent, Santa Clara, CA, USA) and quantified using a fluorometric quantification method (Qubit, USA) that uses dsDNA-binding dyes. Dilutions were done based on the quantified DNA using 10 mM Tris Buffer (pH 8.5). Five microliter (5 µl) aliquots of diluted DNA from each library were mixed for pooling libraries with unique indices. The pooled final DNA library (4 nM) was denatured and sequenced on an Illumina MiSeq System using paired 300-bp reads to generate high-quality reads of the V1–V3 region. Finally, raw fastq files were obtained after trimming the adapters and primer sequences for further bioinformatics analysis.

Data analyses

The obtained raw sequence (Fastq) datasets were initially scrutinized for PCR artefacts and low-quality reads (reads with > 50% bases having a quality score < 2) using ngsShoRT (next-generation sequencing short reads) trimmer as described by Chen et al. (2014). Following the screening process, all the sequence data sets were analysed using the Mothur Pipeline v.1.40.0 as described by Schloss et al. (2009). During the analysis, sequence reads containing low nucleotides (< 50 nts), ambiguities (> 2%), and homopolymers (7%) were excluded, along with sequences of mitochondrial and chloroplast origins. Chimeric sequences were removed using UCHIME algorithm as described by Edgar et al. (2011), while non-chimeric reads were classified using the Naïve Bayesian classifier algorithm as described by Wang et al. (2007) against the SILVA database version 132 (Quast et al. 2013) with a confidence threshold of 80% to assign taxonomic identity of bacteria. A pairwise distance matrix (Euclidean distance matrix) was created from the curated aligned datasets to group sequences into Operational Taxonomic Units (OTUs) at a sequence similarity of 97% for identification. The diversity indices (Shannon–Weaver and Simpson indices) and microbial community richness index (Chao 1) were calculated at the genetic distance of 0.03 to measure the diversity of bacterial species among the data sets. The identified dominant OTUs at phylum and class level were used to generate stacked bar chart using GraphPad prism v 8.01 Software. The raw high-throughput sequencing data were deposited into the NCBI Sequence Read Archive database (SRA accession: PRJNA604090).

Functional prediction analysis

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) software package (Douglas et al. 2019) was used to predict and understand the potential functional capabilities of the bacterial communities at the different sampling sites. Towards this, the Nearest-Sequenced Taxon Index (NSTI) value was used to validate the reliability of predicted gene families and metabolic pathways as NSTI values ≥ 2 are considered noisy and unreliable. PICRUSt2 predicts gene families based on their nearest known taxonomic neighbour as interpolated from available fully sequenced genomes. In this case, the sequences were first aligned to HMMER (Eddy 2011). The aligned sequences were subsequently placed into a reference tree using EPA-NG (Barbera et al. 2019) and gappa (Czech et al. 2019). Normalisation of multiple 16S rRNA gene copies and prediction of gene families were achieved using Castor—a hidden state prediction tool (Louca and Doebeli 2018). The predicted gene families were subsequently collapsed into MetaCyc pathways using MinPath (Ye and Doak 2009) to have a simple overview of the beach sand bacterial community metagenome. Detected pathways were subjected to differential abundance analysis using Aldex2 package (Fernandes et al. 2014) in R statistical software. Pathways with significant p values (≥ 0.05), Benjamini–Hochberg’s FDR score ≤ 0.05, and effect size ≥ 0.6 were considered differentially abundant. The heatmap of the predicted relative abundances of genes related to different functions was generated using heatmap.2 package in R (v3.5.2) (R Core Team 2019).

Results

Heavy metal analysis

Aluminum (Al) was the most abundant metalloid in all sampling sites, with concentrations ranging from 508 ± 54.8 to 2634 ± 187 µg/L (Table 1). Of the heavy metals, Cd was the least abundant metal in all sampling sites (0.1 ± 0 to 0.5 ± 0 µg/L), while Mn was the most abundant (20.9 ± 1.3 to 92.0 ± 4.6 µg/L).

Table 1

Heavy metal concentrations in sand samples from four beach sites along the Indian Ocean coastline in Durban

		Metal concentrations (µg/L)
		B	Al	Cr	Mn	Ni	Cu	Zn	As	Cd	Pb
Limit of detection		0.002	0.045	0.007	0.001	0.015	0.009	0.005	0.053	0.007	0.042
Sampling sites	North beach	5.9 ± 0.3	508 ± 54.8	2.78 ± 0	20.9 ± 1.3	0.8 ± 0.2	0.4 ± 0.6	2.4 ± 0.1	2.9 ± 0.9	0.5 ± 0	0.9 ± 0.2
	South beach	10.2 ± 1.4	2634 ± 187	12.75 ± 0.4	92.0 ± 4.6	5.0 ± 0.1	4.5 ± 0.7	8.9 ± 1.2	6.3 ± 1.2	0.1 ± 0	2.0 ± 0
	Central beach	14.9 ± 2.6	1383 ± 97	11.8 ± 0.9	68.1 ± 2.6	3.4 ± 0.3	2.1 ± 0.4	8.0 ± 0.9	6.1 ± 0.8	0.1 ± 0	2.9 ± 0.1
	Harbour beach	10.5 ± 1.8	1218 ± 84	9.0 ± 0.1	24.8 ± 0.3	2.4 ± 0.4	5.5 ± 0.7	26.7 ± 1.8	2.4 ± 0.4	0.1 ± 0	5 ± 0.4

16S rRNA targeted sequencing analysis

A total of 100 424 raw sequence reads were recovered from all beach sand samples. After cleaning of the sequences, 58 568 sequences were obtained, representing about 58% of the raw sequences. The lowest number of reads (2 215) was obtained in Central Beach sand sample, while the highest number of reads (30 894) were obtained in the North Beach sand sample. A total of 3 204 operational taxonomic units (OTUs) were obtained in all sand samples combined, with 246 OTUs from Central Beach sample and 1 314 OTUs from Harbour Beach sample. The sample-to-sample abundance-based coverage estimator (ACE) and Chao1 indices, which are both used to estimate bacterial species richness, closely matched each other, indicating high levels of accuracy in the analysis. Both Chao1 and ACE indices revealed that the Harbour Beach sample had the highest bacterial species richness followed by the North Beach sample, while the Central Beach sample had the lowest observed species richness. The bacterial species diversity of the sand samples was higher when estimated by the Shannon–Weaver index (H) than when it was estimated using the Simpson’s index (D). The highest bacterial diversity (Shannon index 5.93) was observed in the Harbour Beach sand sample, in which was also recorded the highest number of bacterial phyla and classes (Figs. 2 and 3). Table 2 provides a summary of the bacterial species richness and diversity statistics.

Fig. 2

Phylum-level bacterial community structure in sand samples from four recreational beaches in Durban

Fig. 3

Class-level bacterial community structure in sand samples from four recreational beaches in Durban

Table 2

Summary of the sequence reads, OTUs, and bacterial diversity and richness indices in four beach sand samples

Sample name	Total reads	After filtering	OTUs	ACE	Chao	Shannon	Simpson
Central	3567	2215	246	393.71	423.60	3.52	0.09
Harbour	14,508	7728	1314	1612.93	1503.00	5.93	0.01
North	52,051	30,894	1071	1243.19	1194.20	3.89	0.10
South	30,298	17,731	573	618.11	593.65	4.23	0.04

Analysis of bacterial community distribution at phylum level resulted in the recovery of 11 major bacterial phyla. Sequences belonging to the phylum Proteobacteria were the most dominant in all samples combined, followed by those of the phyla Firmicutes, Actinobacteria, Bacteroidetes, and Gemmatimonadetes. At the Central Beach, 40% of the bacterial sequences obtained belonged to the phylum Firmicutes, 38% belonged to the phylum Actinobacteria, and 20% to the phylum Proteobacteria. The remainder of the recovered sequences belonged to some minor phyla, as represented in Fig. 2. At the Harbour Beach, Proteobacteria was the most dominant phylum with 45% of the recovered sequences. This was followed by the phyla Firmicutes and Actinobacteria each with 14% of the recovered sequences, Bacteroidetes with 10% of the sequences, and Planctomycetes and Verrucomicrobia, each of which was represented by 4% of the recovered sequences. The remainder of the sequences represented some minor phyla represented by less than 1% of the recovered sequences. Samples from the North Beach were dominated by three major phyla, namely, Firmicutes (44% of the recovered sequences), Actinomycetes (14%), and Proteobacteria (39%), while the rest of the sequences represented some minor phyla. Similarly, South Beach sand samples were dominated by sequences representing three major phyla, namely, Proteobacteria (86%), Firmicutes (12%), and Actinobacteria (10%).

At class-level distribution, 16 bacterial classes were identified, the most dominant of which were, in order of overall sequence abundance, Alphaproteobacteria, Bacilli, Actinobacteria, Gammaproteobacteria, Betaproteobacteria, and Clostridia. Bacterial class sequence abundance closely mirrored the phylum-level distribution where most of the dominant classes (≥ 1% sequence representation) belonged to dominant phyla identified in Fig. 2. The Harbour Beach sand samples had the highest bacterial diversity at class level with 14 bacterial classes being represented by sequences with an abundance of ≥ 1% of the total sequences recovered. The most dominant of the classes were in order of magnitude, Gammaproteobacteria, Deltaproteobacteria, Alphaproteobacteria, Clostridia, Acidimicrobiia, Actinobacteria, Flavobacteria, Bacteroidia, Bacilli, Kiritimatiellae, Oligoflexia, Planctomycetia, Phycisphaerae, and Betaproteobacteria (Fig. 3). The Central Beach sand sample harboured the second most diverse bacterial community with 10 bacterial families each having ≥ 1% sequence representation. In order of representative sequence abundance, the identified classes were Actinobacteria, Bacilli, Clostridia, Alphaproteobacteria, Gammaproteobacteria, Betaproteobacteria, Deltaproteobacteria, Gemmatimonadetes (phylum Gemmatimonadetes, Fig. 2), Bacteroidia, and Coriobacteriia. The North Beach sample was dominated by bacteria belonging to the following classes; Bacilli, Alphaproteobacteria, Actinobacteria, Betaproteobacteria, Gammaproteobacteria, Clostridia, Bacteroidia, Physcisphaerae, and Acidimicrobia. The least diverse bacterial community was observed in the South Beach sample, with 5 dominant bacterial classes identified as follows: Alphaproteobacteria, Betaproteobacteria, Bacilli, Gammaproteobacteria, and Actinobacteria.

Genus-level taxonomic analysis showed the presence of 1 163 bacterial genera in all samples. Of these, 32 bacterial genera had a sequence representation of ≥ 1% in at least one sample, while the rest had under 1% sequence representation and were therefore classified as trace members. In terms of ubiquity, bacterial genera belonging to the phylum Proteobacteria dominated, making more than 50% of the recovered genera. However, in terms of percentage sequence abundance per sample, sequences representing the genus Bacillus (phylum Firmicutes) were the most abundant at 37.06% in the North Beach sample, followed by the genus Bifidobacterium (Actinobacteria) at 35.76% and 11.05% in Central Beach and North Beach samples, respectively. The genus Lactobacillus (Firmicutes) was the third most abundant with a sequence coverage of 28.80% in Central Beach sample followed by the genus Methylobacterium (Proteobacteria) at 14.02% sequence coverage in the South Beach sample. The phylum Bacteroidetes was also represented among the top 32 bacterial genera, while sequences representing bacterial genera belonging to other phyla including Planctomycetes, Gemmatimonadetes, Verrucomicrobia, Lentisphaerae, Acidobacteria, Chloroflexi, and Cyanobacteria were present in trace numbers. A summary of the bacterial phyla, class, order, and genera recovered from the beach sand samples is given in Table 3.

Table 3

Bacterial communities (%) at the genus level in recreation sand beach samples

Taxonomy (Phylum; Class; Order; Genus)	Bacterial communities %
	Central	Harbour	North	South
Proteobacteria; Betaproteobacteria; Burkholderiales; Achromobacter	*Tr	Tr	4.43	5.69
Bacteroidetes; Flavobacteria; Flavobacteriales; Actibacter	Tr	1.14	Tr	Tr
Proteobacteria; Gammaproteobacteria; Oceanospirillales; Alkalimarinus	Tr	1.06	Tr	Tr
Firmicutes; Clostridia; Clostridiales; Andreesenia	Tr	1.20	Tr	Tr
Proteobacteria; Betaproteobacteria; Burkholderiales; Aquabacterium	Tr	Tr	2.36	3.63
Actinobacteria; Actinobacteria_c; Micrococcales; Arthrobacter	Tr	Tr	Tr	1.32
Proteobacteria; Alphaproteobacteria; Rhodospirillales; Azospirillum	Tr	Tr	Tr	1.21
Firmicutes; Bacilli; Bacillales; Bacillus	Tr	Tr	37.06	2.54
Actinobacteria; Actinobacteria_c; Bifidobacteriales; Bifidobacterium	35.76	2.78	11.05	2.43
Proteobacteria; Alphaproteobacteria; Rhizobiales; Bradyrhizobium	Tr	Tr	2.33	5.23
Firmicutes; Clostridia; Clostridiales; Clostridium	4.42	1.77	Tr	Tr
Proteobacteria; Deltaproteobacteria; Desulfobulbaceae_o; Desulfofustis	Tr	3.96	Tr	Tr
Proteobacteria; Deltaproteobacteria; Desulfobacterales; Desulfosarcina	Tr	2.14	Tr	Tr
Proteobacteria; Gammaproteobacteria; Enterobacterales; Enterobacteriaceae_g	1.22	Tr	4.97	6.27
Firmicutes; Clostridia; Clostridiales; Gottschalkiaceae_uc	Tr	2.25	Tr	Tr
Firmicutes; Bacilli; Lactobacillales; Lactobacillus	28.80	1.73	4.27	1.21
Proteobacteria; Alphaproteobacteria; Rhizobiales; Methylobacterium	Tr	Tr	4.64	14.02
Actinobacteria; Actinobacteria_c; Frankiales; Nakamurella	Tr	1.40	Tr	Tr
Actinobacteria; Actinobacteria_c; Propionibacteriales; Nocardioides	Tr	Tr	Tr	3.40
Proteobacteria; Betaproteobacteria; Burkholderiales; Noviherbaspirillum	Tr	Tr	Tr	3.20
Proteobacteria; Gammaproteobacteria; Oceanospirillales; Oceanospirillum	Tr	2.23	Tr	Tr
Proteobacteria; Alphaproteobacteria; Rhizobiales; Ochrobactrum	1.17	Tr	1.82	2.25
Proteobacteria; Alphaproteobacteria; Rhodobacterales; Phaeobacter	1.13	Tr	Tr	Tr
Proteobacteria; Betaproteobacteria; Burkholderiales; Ralstonia	Tr	Tr	1.06	3.07
Proteobacteria; Alphaproteobacteria; Rhizobiales; Rhizobiaceae_uc	Tr	Tr	Tr	1.21
Proteobacteria; Alphaproteobacteria; Rhizobiales; Rhizobium	Tr	Tr	Tr	1.24
Proteobacteria; Alphaproteobacteria; Rhodobacterales; Rhodobacteraceae_uc	Tr	Tr	Tr	8.18
Bacteroidetes; Flavobacteria; Flavobacteriales; Robiginitalea	Tr	1.09	Tr	Tr
Proteobacteria; Alphaproteobacteria; Sphingomonadales; Sphingomonas	Tr	Tr	Tr	1.75
Firmicutes; Clostridia; Clostridiales; Sporobacter	1.31	Tr	Tr	Tr
Proteobacteria; Alphaproteobacteria; Rhodospirillales; Terasakiella	Tr	7.76	Tr	Tr
Proteobacteria; Gammaproteobacteria; Enterobacterales; Wigglesworthia	1.35	Tr	Tr	Tr
Others	17.88	66.65	20.38	31.52

*Tr trace, representing sequence coverage of less than 1%

Canonical correspondence analysis (CCA)

Determination of the effects of heavy metal concentration on distribution and prevalence of bacterial metacommunities by CCA has showed that bacteria were differentially distributed with respect to heavy metal concentrations (Fig. 4). At the harbour sampling site, for instance, the metals Pb and Zn influenced the distribution of Flavobacteria, Bacteroidia, and Deltaproteobacteria, while the distributions of Clostridia and Gammaproteobacteria were correlated to the concentrations of barium (B) and Cr, respectively. Bacterial distribution at the South and North beach sampling sites, largely Alphaproteobacteria and Betaproteobacteria, was correlated to the metals Ni, Mn, As, and Al. However, bacterial distribution at Central sampling site was not statistically linked to the concentration of any metal species.

Fig. 4

Canonical correspondence analysis (CCA) showing the effects of different heavy metal concentrations on bacterial distribution patterns along the four sample sites

PICRUSt2 prediction of functional genes among the sand-borne bacterial metacommunities revealed the presence of both pathogenic and antibiotic resistance pathways. Identified antibiotic resistance genes included the peptidoglycan biosynthesis II (staphylococci) exhibited in the Harbour, North, and South samples (Fig. 5). Others included the peptidoglycan biosynthesis III (mycobacteria), IV (Enterococcus faecium), and V (beta-lactam resistance), as well as the polymyxin resistance gene, which was identified in all samples. The identified virulence genes included the sitA gene, the fimB gene, the aerobactin synthase gene, and pilL gene, all of which were detected in bacterial metacommunities from all sampling sites.

Fig. 5

Virulence enzymes and pathways detected in the beach samples (differences in the overall abundance of virulence and resistance genes in the various sample were significant; p = 0.028)

Discussion

Data in Table 2 show that there were slight differences between the ACE and Chao1 species richness estimators. Hughes et al. (2001) point out that the ACE and Chao1 estimators are similar in that they are nonparametric. They are different in that the Chao1 estimator gives greater weight to low abundance species (species with less than 10 individuals in a sample), while the ACE estimator gives greater weight to species sample coverage (species with more than 10 individuals in a sample) (Kim et al. 2017). In this study, however, both estimators still showed high bacterial species richness in all beach sand samples, as also observed in other beach sand samples (Mudryk et al. 2013; Romão et al. 2017). This confirms earlier findings that beach sands are microhabitats teaming with microbial life (Sabino et al. 2014; Whitman et al. 2014; Solo-Gabriele et al. 2015). The disparities observed between the Shannon–Weaver diversity and the Simpson’s diversity indices imply that, while beach sand samples had high species richness, they had low species evenness. This trend has also been observed in other studies of a similar nature, including that of Gobet et al. (2012).

While most microbes found in beach sand are harmless, some are linked with human disease (Rettner 2018). Therefore, understanding the microbial community structure and biodiversity of beach sand is useful for assessing the ecological health of the beach sand, as higher biodiversity usually confers better resistance and resilience against environmental perturbation and contamination (Cui et al. 2013). The dominance of bacteria belonging to the phyla Proteobacteria, Bacteroidetes, Planctomycetes, Actinobacteria, and Gemmatimonadetes in beach sands has earlier been documented in coastal beach samples (Zheng et al. 2014; Romão et al. 2017). This may indicate that these phyla have evolved to be the main indigenous bacterial communities in coastal intertidal zones. Bacteria belonging to the phylum Proteobacteria are known for their metabolic plasticity (Esposti 2014), causing them to be able to colonise the most diverse environments compared to other bacterial phyla. This could be the reason why, at genus level, greater than 50% of recovered OTUs belonged to the phylum Proteobacteria. As well as being metabolically diverse, this phylum is known for containing the most number of bacterial pathogens in its ranks (Rizzatti et al. 2017).

However, many other bacterial phyla identified in our study house essentially non-pathogenic bacteria. Most such bacteria have fundamental roles in the biogeochemical cycles of nutrients/mineral in the sediments (Zheng et al. 2014; Cardenas 2016). For example, bacteria belonging to the phylum Chloroflexi have been recognized for their key role in oil degradation (Zheng et al. 2014), while bacteria of the genus Ochrobactrum are known for hydrocarbon degradation (Octaviany et al. 2019).

Sequences representing harmless bacterial populations were the most abundant in this study. These included sequences for Bifidobacterium longum, Lactobacillus paracasei, Lactobacillus helveticus, Bifidobacterium animalis, Bifidobacterium pseudolongum, Bifidobacterium bifidum, and Lactobacillus buchneri, which are largely used in the food industry as probiotics (Margolles and Sánchez 2012; Smokvina et al. 2013; Sugahara et al. 2015; Wong et al. 2019). Nonetheless, because these probiotic bacteria are natural flora of the gastrointestinal tracts of humans, and that of other animals (Mikkelsen et al. 2003), the abundance of sequences representing these bacterial species in beach sand samples serves as an important indicator of faecal pollution. Such faecal pollution could emanate either from faecal shedding by sand bathers, direct faecal deposition by animals and birds, or from faecally contaminated sea water washing onto the beach sands (Whitman et al. 2014). However, the presence of other Bifidobacterium sp. like Bifidobacterium dentium can be used to determine faecal contamination of exclusively human origin (Nebra et al. 2003; Furet et al. 2009), and can be useful in determining the public health safety of recreational beaches. Although such allochthonous microorganisms get subjected to various biotic and abiotic pressures that affect their fate in the new environment (Feng et al. 2010; Gobet et al. 2012), microorganisms of enteric origin in beach sands may be incidentally ingested, leading to carriage, colonization, or even infection.

Sequences belonging to the genera Achromobacter xylosoxidans and Methylobacterium were both detected in the South and North beach sand samples. Achromobacter xylosoxidans is an emerging, multidrug-resistant opportunistic pathogen which invades cystic fibrosis (CF) patients’ airways, as well as cause a wide variety of infections in immunocompromised patients (Marion-sanchez et al. 2019). Its presence in recreational beaches presents a potential health risk to beach goers, more so to immunocompromised persons. Amoureux et al. (2013) earlier reported that this bacterium is innately resistant to cephalothin, cefoxitin, cefotaxime, aztreonam, and aminoglycosides, and frequently shows acquired resistance to carbapenems, ceftazidime, and ciprofloxacin, drastically limiting therapeutic choices for infected individuals. However, that largely depends on the environment from which it is isolated, with clinical strains having proven pathogenesis. Apart from this study, there are no previous reports linking A. xylosoxidans to beach sands. Previous reports show that it has been found in some plants, polluted soils, well water, domestic and hospital drains, as well as freshwater bodies frequently used for recreational purposes (Amoureux et al. 2013). Being an emerging pathogen, the natural habitat for A. xylosoxidans, as well as its medium of spread are not yet known. Therefore, the identification of potential environmental reservoirs of this bacterium might aid in the prevention of infection among CF patients as well as immunocompromised individuals. Transmission of infectious diseases in beach environments can occur via direct exposure to microbes found in the sand through such routes as dermal contact, contact with eyes and ears, inhalation, and ingestion (Solo-Gabriele et al. 2015). The varied, though inevitable interactions between different beach sand zones and beach water make beach sand a potential source of such pathogens in beach water also (Feng et al. 2010; Cui et al. 2013), which increases the risk of exposure to not only the sand bathers but the surfers and swimmers too. Methylobacterium sp., meanwhile, is an opportunistic bacterial pathogen in immunocompromised persons. It is reported to form biofilms, and to be tolerant to disinfecting agents, high temperatures, and drying (Kovaleva et al. 2014), a befitting reason as to why their sequences were found in relatively high percentage abundances in beach sand samples in this study.

We also recovered sequences belonging to the species Aeromicrobium erythreum from all sampling sites though in low relative abundances ranging from 0.06% of the total sequences recovered from the Harbour Beach sample to 0.41% of those recovered from the Central Beach sample. Bacteria belonging to the genus Aeromicrobium have previously been isolated from sea water by Bruns et al (2003) who described it as an obligately salt-dependent Gram-positive bacterium affiliated to the family Nocardioidaceae, within the order Actinomycetales. While there are currently no reports of its recovery from other environments except sea water, the recovery of sequences representing a species in this genus from beach sand provides further proof of tide-assisted microbial interactions between the beach sand and the sea water. Alternatively, this finding could suggest that this bacterium could be found in more diverse environments, and not just in the seawater. However, both the health significance and the influence of anthropogenic pressure on the distribution of this bacterium is currently unknown. Further studies are needed to, at the least, establish its ecological significance.

PICRUSt2 determination of functional genes also showed that beach sand bacterial metacommunities are important reservoirs of antibiotic and virulence genes. And, just like in other studies (Sandaa et al. 1999; Yao et al. 2016), canonical correspondence analysis (CCA) in this study revealed differential but positive correlations between heavy metal concentrations and certain bacterial metacommunities. The relationship between antimicrobial resistance and potentially toxic metal resistance in bacteria has been a subject of study for a long time (Chen et al. 2019; Quero et al. 2015). Findings suggest that environmental heavy metal concentrations both shape the microbial community compositions as well as induce antimicrobial resistance by either cross- or co-resistance phenomena (Chen et al. 2015; Nguyen et al. 2019; Yazdankhah et al. 2018). Environmental reservoirs of antibiotic resistance genes such as those observed in this study are of particular public health concern, considering that polymyxin is an antibiotic of last resort in the treatment of Gram-negative bacterial infections (Srinivas and Rivard 2017; Li et al. 2019). By almost the same measure, β-lactams are the first and most frequently used class of antibacterial agents used to treat severe infections due to Gram-positive bacteria (Bugg and Walsh 1992; Mainardi et al. 2005). However, peptidoglycan biosynthesis II in staphylococci and peptidoglycan biosynthesis IV in enterococci usually give rise to β-lactams resistance due to production of low-affinity penicillin-binding proteins (PBPs), like the PBP2a D,D-transpeptidase protein which results in methicillin resistance in staphylococci, and the PBP5 D,D-transpeptidase protein which results in ampicillin resistance in enterococci (Mainardi et al. 2000, 2008). Furthermore, some enterococci strains have been found to acquire antibiotic resistance by completely bypassing the D,D-transpeptidase by an L,D-transpeptidase which confers resistance towards most β-lactam antibiotics on the organism (Mainardi et al. 2008).

Type IV Pili, that were possessed by the bacterial metacommunities from all sampling sites, are important virulence factors used by many pathogens including Pseudomonas aeruginosa for attachment to surfaces and host tissues and twitching motility (Kilmury and Burrows 2016). Transcriptional activation of the fimB gene in uropathogenic Escherichia coli (UPEC) leads to increased expression of type 1 pili, a chief virulence factor responsible for UPEC pathogenicity (Zhang et al. 2016). Urinary tract infections afflict mostly women, with UPEC emerging as the primarily agents of infections in humans (Schwan et al. 2018). Inside their hosts, pathogenic bacteria are often have to contend with extremely low bioavailability of iron, which thus becomes a limiting factor for survival (Bailey et al. 2018). Bacterial pathogens circumvent this limitation using unique strategies to scavenge iron, including the synthesis, secretion, and reuptake of iron chelators (siderophores) such as aerobactin, which has since been demonstrated to be critical for virulence in pathogens like Klebsiella pneumoniae (Bailey et al. 2016) and Vibrio mimicus (Moon et al. 2004). The sitA gene is most prominent among bacterial pathogens which utilize it to mobilize iron and manganese inside eukaryotic cells and is therefore associated with bacterial virulence (Runyen-Janecky et al. 2003; Tivendale et al. 2009). Environmental bacteria harbouring reservoirs of both virulence and antibiotic resistance genes, therefore, pose increased risks of human infections, particularly in settings like beaches, since beach visitors tend to spend more time in contact with beach sand than with water (Whitman and Nevers 2003). Chances of infection are significantly increased if higher pathogen densities are found in the sand than in the water column, as previously observed (Sato et al. 2005). In particular, children are exposed to greater risk, because they spend more time playing and digging in the sand, where potentially pathogenic bacteria are likely to persist longer due to adsorption to sand particles, unlike free bacteria in the water. In what is a limitation of our study, PCR confirmation of the presence of virulence and antibiotic resistance genes as predicted by PICRUSt2 results could have made our findings very concrete.

The identification of bacterial sequences belonging to the family Enterobacteriaceae in all but the Harbour Beach sample could indicate the likely contamination of beach sands with enteric microorganisms, which might translate to increased health risks to beach goers as most of these microbes are potential pathogens. Moreover, most environmental pathogenic strains of this family are known to be multi-drug-resistant (Gonçalves et al. 2019).

Sequences belonging to some unique bacterial genera were also found, some of them with no previous history of being isolated from beach sands or marine environments. These included sequences belonging to the genera Nakamurella, Robinginitalea and Teresakiella (Harbour Beach sample), Phaeobacter, Ralstonia, Sporobacter, and Wigglesworthia (Central Beach sample). Bacteria of the genus Wigglesworthia are otherwise known to reside within the cytoplasm of differentiated epithelial cells (bacteriocytes) of tsetse flies (Pais et al. 2008; Soumana et al. 2014). Bacteria belonging to the genus Ochrobactrum and the family Rhodobacteraceae have been identified as etiological agents of the black band disease in corals (Miller and Richardson 2010), while Teresakiella sp. (either pusilla or brassicae) is thought to colonise marine shellfish (Han et al. 2016), and could have been deposited onto the beach sand by ocean tides.

In conclusion, the results of this study demonstrate that beach sand-borne bacteria are potential reservoirs of virulence and antibiotic resistance genes. Contamination of beach sands with heavy metals selects for both heavy metal resistance and antibiotic resistance in beach sand microbial communities. The abundance of pathogenic bacterial OTUs in beach sand shows the likelihood of health risks to people engaging in recreational activities in beach sands. In particular, immunocompromised persons are at an increased risk of contracting bacterial infections after beach visits, especially taking into consideration the likelihood of incidental ingestion of beach sand and water. To validate this possibility, future studies may need to quantitatively assess the microbial risk that beach goers are exposed to. However, not all microorganisms in beach sand are pathogenic as others play significant ecological roles including suppression of bacterial pathogens, nutrient recycling, and also as bio-filters.