Freshwater sponges in the southeastern U.S. harbor unique microbiomes that are influenced by host and environmental factors

Site description, sample collection, and nutrient analysis

This dataset began with pilot 16S rRNA gene cloning data from the first year of sponge collections (2016) in this region (Supplemental Information), which motivated a larger geographic search for sponges the following year (2017) and a shift to high throughput sequencing. Because of this, three additional sampling sites were added throughout the year yielding an uneven sampling design and different timing of collection at sites. Sponges and corresponding water samples were collected from a total of four locations near Boone, North Carolina, USA in 2017 (Fig. 1; Table 1): New River at the New River State Park in July, August, and twice in September (NR; collection permit 2021_0548), Watauga River in September (WR), Jacob Fork River in November (JF), and a mussel hatchery pond and indoor tank with mussels in September (Hatch). The hatchery is at a similar elevation to the Jacob Fork River site and while not directly connected, they are in the same water basin (Catawba River basin). In contrast, the New River and Watauga River sites are each in different water basins (New and Watauga River basins, respectively).

All sites were surrounded by a mix of agriculture and woodlands and generally close to a road, which means that all sites were subject to sporadic road runoff as well as agricultural runoff. The largest physical or geographic difference between the sites is elevation as described above and the corresponding air temperature which tends to be several degrees warmer at the lower elevation sites. Sponges in the New, Watauga, and Jacob Fork rivers were found on the underside or sides of rocks that were large enough to not tumble easily with the current. Sponges were never found on small rocks, pebbles, or wood debris at these sites. The rivers of New, Watauga, and Jacob Fork are all fourth order streams and stream flow rate and depth vary considerably over time. At the times of sampling, the water depth generally varied between 0.25 and 1 m depth. After rain events the water frequently rises to 1.5–2 m depth at these sites and the flow rate is generally increased, while after periods of drought the water may be barely covering the rock substrate. The mussel hatchery pond is seasonal and spring fed from the surrounding Catawba River watershed. The sponges collected from the pond were from approximately 0.75 m depth growing on bricks used to weight the mussel aquaculture containers. The catchment areas, calculated by the United States Geological Survey (streamstats.usgs.gov/ss/) for each sampling site are: 606.1 km² (New River), 244.8 km² (Watauga), 213.4 km² (Jacob Fork), and 1.4 km² (mussel hatchery).

To address the hypotheses that sponge microbiomes would differ 1) from the surrounding water, and 2) across study locations, sponge and water samples were collected at the same time each time of sampling (except for early-September at the New River; Table 1). Individual sponges were photographed, then scraped with a razor blade and subsections of the sponge were preserved in 80% ethanol and in DNA buffer (Seutin, White & Boag, 1991) and stored on ice until returning to the laboratory. Water samples (500 mL each, n = 2–3 per sampling point) were collected in acid-cleaned glass bottles and each bottle was rinsed in the river water three times before collecting water for analysis. The water was collected near the surface in arbitrary locations across the width of the river. The mussel hatchery site contained a small pond that was used to grow mussels on nets and an indoor flow through system with pond water to maintain mussels in tanks: five sponge specimens were collected from the pond and three from the wall of indoor tanks, alongside one tank water and two pond water samples. The hypothesis that sponge and water microbiomes would differ by site was also addressed by this sampling scheme of paired sponge and water samples at each of the four sites.

To address the hypothesis of increased similarity in the sponge and water microbiomes during the growth of the sponge, at one site, the New River, we sampled four times over the growing season (July–September) and we attempted track individual sponges over time. Attempts were made to mark the rocks that sponges were collected from using flagging tape, nail polish, and spray paint, but none of these lasted in the water. Thus, as best that we can tell (based on photos of the sponge and rock), the only resampling of the same sponge individuals were the final three individuals collected from the New River in middle September, NR41, NR42, NR43, that are the same individuals as NR28, NR30, and NR29 (sampled in early September), respectively.

In the laboratory, the ethanol in the sample tubes was replaced with fresh 80% ethanol and stored at 4 °C and sponges in DNA buffer were stored at −20 °C until DNA extraction. Water samples were filtered through a 0.22 µm polyethersulfone (PES) filter (Whatman, USA), the filter was preserved in DNA buffer and stored at −20 °C and the filtrate was frozen until used for inorganic nutrient analysis. Additional water samples were collected in 2018 for inorganic nutrient analysis and were filtered in the same process. The filtrate from both years was used for inorganic nutrient analysis. The inorganic ions, nitrate, phosphate, sulfate, and chloride were measured using a Thermo Scientific (formerly Dionex) ion chromatograph in the Chemistry and Fermentation Sciences department at ASU equipped with an IonPac AS11-HC column set (analytical column 4 × 250 mm; guard column 4 × 50 mm) and a conductivity detector. The background conductivity was suppressed using an AERS 500 (4 mm) suppressor operated at 112 mA. The eluent was 25 mM KOH at a flow rate of 1.4 mL min⁻¹ under isocratic conditions. The injection volume was 20-μL and the column temp 30 °C. A standard mix at 1 ppm (Environmental Express, Charleston, SC) was used to create a standard curve for quantification of each compound. The water samples were analyzed in duplicate (technical replicates) to obtain multiple measurements. The values from the technical replicates were averaged and then sample replicates were averaged and standard deviation calculated for visualization.

Water temperature and pH were not recorded at the time of sampling in 2017. However, these data are publicly available for Cranberry Creek (36.5697, −81.327), a tributary of the New River, through the New River Conservancy (https://newriverconservancy.org/) and were used here as an indicator for changes in temperature and pH with season as opposed to absolute values of temperature and pH.

Sponge species identification was accomplished by spicule analysis. Spicule preparations were conducted by dissolving a small piece of sponge in bleach overnight in a microfuge tube. The spicules were rinsed with distilled water, then mounted on a slide to view with a compound microscope. The Ecology and Classification of North American Freshwater Invertebrates (Reiswig, Frost & Ricciardi, 2010) was used as a key for identifying the sponges. Megascleres, microscleres and gemmuloscleres (when present), were observed and used to identify the sponge.

DNA extraction and 16S rRNA gene sequencing

DNA extractions from subsections of sponge and the whole filter were performed using the cetyltrimethylammonium bromide (CTAB) method originally designed for plant DNA extractions (Lipp et al., 1999). Extracted DNA was checked for concentration and quality using a Nanodrop spectrophotometer ND1000 and sent to Nova Southeastern University for sequencing on an Illumina MiSeq. Polymerase chain reaction (PCR) was performed following published protocols from the Earth Microbiome Project (http://www.earthmicrobiome.org/protocols-and-standards/). PCR reactions were run using barcoded 16S rRNA primers (515 F and 806 R: (Caporaso et al., 2011; Walters et al., 2015)). Barcoded PCR products were cleaned using AMPure beads (Beckman Coulter), and concentrations of the cleaned products were measured using a Qubit fluorometer. Cleaned products were diluted to 4 nM and then pooled in equal volumes. Sample preparation and loading onto the Illumina MiSeq followed a standard Illumina protocol, with the exception of using custom sequencing primers from the EMP protocol. Samples were sequenced using a 500 cycle V2 chemistry kit (Illumina) for 16S rRNA gene sequencing to obtain 250 bp single-end amplicons. In additional to 16S rRNA gene analysis, two marker genes of the sponge (cytochrome oxidase I (COI) and 18S rRNA) were amplified and archived (Supplemental Information), but due to poor representation of the freshwater sponge genes in the databases, these sequences were not useful to identify the sponges in this study.

Data analysis

The sampling map was created using ESRI ArcMap v10.8.1 Online (2024). Raw FastQ files were analyzed with the Dada2 pipeline (Callahan et al., 2016) within the R statistical program (R Core Team, 2021). Dada2 was used to trim, filter, and deduce representative amplicon sequence variants (ASVs). Within the Dada2 pipeline, reads were trimmed to 150 (forward) and 180 (reverse), with average quality scores between 20 and 30 for the remaining reads and were denoised and filtered for chimeric sequences and singletons (one read for an ASV) in Dada2. The SILVA database v138 was used as a reference ribosomal gene database. The Silva SSU database was used for taxonomic classification (v138; Quast et al., 2013). Phyloseq (v1.36.0; McMurdie & Holmes, 2013) and ggplot2 (Wickham, 2016) were used for downstream analysis and visualization.

Alpha and beta diversity analysis was performed with the vegan package (Oksanen et al., 2022) in R. Alpha diversity was compared across sample type and location with analysis of variance (ANOVA) following tests for assumptions of ANOVA and followed by Tukey’s Honestly Significant difference post-hoc tests for all locations except for Watauga river sponges. Beta diversity analyses used relative proportions of reads for non-metric multidimensional scaling (nMDS, Phyloseq ‘ordinate’ function) and comparison of taxa across samples. Ordination by nMDS used Bray-Curtis or weighted Unifrac (Lozupone & Knight, 2005) distance matrices. Differences in composition of taxa across sample groups (e.g., sponge vs. water samples) was deduced using a permutational analysis of variance (PERMANOVA) with the Adonis2 function in the vegan package in R. Identification of ASVs of interest between sponge and water samples (i.e., taxa present or present in high proportion in one but not the other sample type) were discovered by (1) heatmap visualization of taxa abundance in each sample type, and (2) the vegan similarity percentage analysis (SIMPER) function in R comparing sample type.

An overview of the sponge and water microbiomes was conducted for all samples (i.e., taxonomic profile and nMDS visualization) and then a more in-depth analysis across month was performed for the New River samples. These included visualization by ordination using nMDS (using Bray-Curtis and Unifrac distance matrices and the Phyloseq ‘ordinate’ function), heatmap visualization (‘plot_heatmap’ function from Phyloseq). PERMANOVA (Adonis2 function) was used to test for differences in the composition of symbiont taxa by month.

Differences in temperature over time at the New River site were assessed with a Kruskal-Wallis rank sum test in R (Hollander & Wolfe, 1973). Differences in inorganic nutrient concentrations between sites or over time were not conducted due to limited replicates.

Data availability

The sequence data are available at NCBI in BioProject PRJNA988097 (16S rRNA gene microbiome sequences), BioProject PRJNA975882 (sponge target loci sequences), and accession PP930792 for the 16S rRNA gene clone. The sequences are also available as Supplemental Files (sponge CO1, 18S rRNA and ITS, and one 16S rRNA gene sequence clone) and at the OSF project site for the ASV sequences (https://osf.io/rpbz7/). The raw environmental data as well as the R script and output files from Dada2 processing are available online through OSF (https://osf.io/rpbz7/).

Sponge identification

Morphological identification using spicules indicated a total of three species. All sponges from the Watauga River and the New River, except for one individual in the New River, were identified as Radiospongilla crateriformis (Potts, 1882). The megascleres were oxeas of one size class with subtle spines (Figs. S1A–S1C, S1G). The gemmuloscleres, when present (Table 1), were birotula with curved spines on each rotule and several spines at the terminal ends of the rotule shaft (Figs. S1A, S1B). For individuals without gemmuloscleres (Table 1), the high similarity in megascleres suggested that these are likely the same species. One individual of the genus Trochospongilla (Fig. S1I) was also observed at the New River site, with the distinctive spiny oxeas, identical to the oxea megascleres of Trochospongilla horrida (Weltner, 1893) identified from the Jacob Fork River (Fig. S1D). The Trochospongilla from the New River did not have gemmuloscleres and thus we refer to it as Trochospongilla sp., although it is likely T. horrida as we observed only T. horrida in the region. The distinctive gemmuloscleres of T. horrida from Jacob Fork were smooth and similar sized birotules (Fig. S1D). Trochospongilla horrida at Jacob Fork was present as small smooth individuals (not shown) and with a convoluted phenotype (Fig. S1L) and both tan color and green suggestive of algal symbionts. The only species observed at Jacob Fork was T. horrida. At the mussel hatchery, the only species observed was Eunapius fragilis (Leidy, 1851) (Figs. S1K, S1J). The hatchery individuals all had smooth and stout oxeas as megascleres and two spiny strongyle gemmuloscleres (Figs. S1E, S1F). The COI genes were nearly identical across all individuals and both the COI and 18S rRNA sequences yielded equivocal hits to multiple other species in the NCBI database (Tables S1, S2). Thus, the genetic information was not used any further and we relied on the spicule-based identifications.

Microbial community composition in sponges and river water

For high throughput amplicon sequencing, a total of 3,287,273 reads (average of 67,087 reads per sample) were obtained from all samples. Filtering low quality and chimeric reads resulted in 1,262,913 total reads (average of 25,773 per sample) that comprised 1,359 amplicon sequence variants (ASVs).

The alpha diversity of sponge samples was higher on average than for water samples (Fig. S2, ANOVA, F_7,41 = 16.84, p < 0.001, Shannon diversity; F_7,41 = 12.32, p < 0.001, Inverse Simpson diversity). Post-hoc tests did not yield significant differences in either diversity metric for sponge vs. water samples at a given site (Tukey HSD, adjusted p > 0.05). However, post-hoc tests (Tukey HSD) did yield significant differences in both diversity metrics between sites for sponge samples, with the hatchery sponges having higher diversity than Jacob Fork (adjusted p < 0.001 Shannon and Inverse Simpson) and New River (adjusted p < 0.001 Shannon and Inverse Simpson), New River sponges were higher in diversity compared to Jacob Fork sponges for Shannon diversity but not for Inverse Simpson (adjusted p = 0.02 Shannon, p = 0.26 Inverse Simpson). Alpha diversity of water samples followed the same pattern as the hatchery water samples, having a higher diversity compared to Jacob Fork (adjusted p < 0.001) and New River (adjusted p < 0.001). The New River water samples had a higher Shannon but not Inverse Simpson in diversity compared to Jacob Fork water samples (adjusted p = 0.02 Shannon, p = 0.87 Inverse Simpson).

For ASVs over 0.01% abundance, these were classified to 13 bacterial phyla. The phyla were dominated by Actinobacteria, Bacteroidota, Verrucomicrobia, Pseudomonodota, and to a lesser extent, Spirochaeta, Patescibacteria, and Chloroflexi. At the class level of taxonomy, there were also only a few taxonomic groups that comprised the most abundant ASVs particularly in the sponge samples (Fig. 2). The predominant classes in sponge samples relative to water included Verrucomicrobiae (phylum Verrucomicrobia), Parcubacteria (Patescibacteria), Leptospirae (Spirochaetota) for the Jacob Fork location, Kapabacteria (Bacteroidetes/Chlorobi group) for the hatchery location, and Alphaproteobacteria (Pseudomonodota) and unclassified ASVs across sponge samples (Fig. 2). There was a higher proportion of Actinobacteria (Actinobacteriota) and Bacteroidota (Fig. 2) in water sample communities, and similar relative abundance of Gammaproteobacteria (Pseudomonodota) and most of the low abundance classes (Fig. 2).

Differences in microbiome composition across sample type and site

Ordination by both ecological (Bray-Curtis) and phylogenetic (weighted Unifrac) distances yielded separation in nMDS space between sponge and water microbiomes among sites (Figs. 3A, 3B). The difference in sample type (sponge vs. water) was supported by PERMANOVA (F_1,48 = 3.1, R² = 0.06, p = 0.003, Bray-Curtis; F_1,48 = 4.38, R² = 0.08, p = 0.003, Unifrac). There was also a significant effect of site for sponges (PERMANOVA, F_3,32 = 15.7, R² = 0.62, p = 0.001,) and for water samples (PERMANOVA, F_3,15 = 16.2, R² = 0.8, p = 0.001) (Figs. 3A, 3B). We also note that the sites tended to be species-specific and were sampled at different times of year, preventing us from testing for an effect of species or time in the microbiome composition across all of our sites. SIMPER analysis yielded 360 ASVs that explained ~80% of variation between sponge and water samples. Of the 360 ASVs, 17 ASVs explained ~30% of variation between sponge and water samples and these 17 were examined further for their taxonomic identification (four taxonomic classes: Bacteroidia, Actinobacteria, Gammaproteobacteria, Verrucomicrobiae, Table 2).

Seasonal variability of the microbiome at the New River site

Additional nMDS ordinations examining New River samples over the months of July, August, and September also support separation by sample type and by month based on Bray-Curtis and Unifrac distance matrices (Figs. 4A, 4B). PERMANOVA supported a significant effect of sample type for New River samples (PERMANOVA, F_1,23 = 4.75, R² = 0.18, p = 0.002 (Bray-Curtis), F_1,23 = 10.9, R² = 0.33, p = 0.001 (Unifrac distance)) (Figs. 4A, 4B). PERMANOVA also supported a significant effect of month for New River samples (PERMANOVA, F_2,23 = 6.8, R² = 0.39, p = 0.001 (Bray-Curtis), F_2,23 = 3.5, R² = 0.25, p = 0.021 (Unifrac distance). A group of three sponge samples collected in late September were unique from the early September samples in that they contained gemmules and the sponge tissue itself was degrading at the time of collection; these three samples grouped together in both New River ordinations (Figs. 4A, 4B). Because Radiospongilla crateriformis accounted for all individuals except for one Trochospongilla horrida individual at the New River site, we could not compare species within this site. It is notable, however, that the T. horrida sample was separated from the R. crateriformis individuals on the ordinations (Figs. 4A, 4B).

At the ASV level for the New River samples, the most abundant ASVs (i.e., those with the most reads), included taxa that were present in both sponge and water and some taxa present only in sponges (Fig. 5 “top 50” ASVs, Fig. S4 “top 100” ASVs). More specifically, a group of taxa present in sponges and not the water samples were only present in the August and September samples, not the July samples (Figs. 5, S4, Table 3). The groups Alphaproteobacteria and Bacteroidetes accounted for the majority of these sponge-derived taxa, with a few Gammaproteobacteria and Actinobacteria (Table 3). BLAST (Altschul et al., 1990) analysis with the sponge-derived Dada2 ASVs yielded similar identity (using the “top” BLAST match) to uncultured freshwater Bacteroidetes and Alphaproteobacteria groups, including bacteria from other freshwater sponges Spongilla lacustris, Lubomirskia baicalensis, and Swartschewskia papyracea, the latter two of which are endemic to Lake Baikal (Table 3).

Environmental variables across sites or season

The pH of Cranberry Creek (near the New River) was stable during 2017 at 6.7 for July through October and over the years from 2013 to 2021 at ~6.5–6.7. Temperature data for Cranberry Creek did not vary significantly by month from 2013–2017 within the months of collection (July, August, September) (Kruskal-Wallis, chi-squared = 3.0526, df = 2, p = 0.2173).

Inorganic nutrients were variable but generally low for chloride, nitrate, and sulfate at New River and Jacob Fork sites (Fig. S5). Phosphate was below detection in all samples. While limited replicates precluded statistical analysis, this preliminary dataset does not suggest any major differences in chloride, nitrate, and sulfate between the New and Jacob Fork sites or over time at the New River site (Fig. S5).