Bloom-forming gelatinous zooplankton occur circumglobally and significantly influence the structure of pelagic marine food webs and biogeochemical cycling through interactions with microbial communities. During bloom conditions especially, gelatinous zooplankton are keystone taxa that help determine the fate of primary production, nutrient remineralization, and carbon export. Using the pelagic tunicate Dolioletta gegenbauri as a model system for gelatinous zooplankton, we carried out a laboratory-based feeding experiment to investigate the potential ecosystem impacts of doliolid gut microbiomes and microbial communities associated with doliolid fecal pellets and the surrounding seawater. Targeted metabarcoding (16S rRNA genes recovering Bacteria/Archaea) and qPCR approaches were used to characterize microbiome assemblages. Comparison between sample types revealed distinct patterns in microbial diversity and biomass that were replicable across experiments. These observations support the hypothesis that through their presence and trophic activity, doliolids influence the structure of pelagic food webs and biogeochemical cycling in subtropical continental shelf systems where tunicate blooms are common. Bacteria associated with starved doliolids (representative of the resident gut microbiome) possessed distinct low-biomass and low-diversity microbial assemblages, suggesting that the doliolid microbiome is optimized to support a detrital trophic mode. Bacterial genera Pseudoalteromomas and Shimia were the most abundant potential core microbiome taxa, similar to patterns observed in other marine invertebrates. Exploratory bioinformatic analyses of predicted functional genes suggest that doliolids, via their interactions with bacterial communities, may affect important biogeochemical processes including nitrogen, sulfur, and organic matter cycling.

Collection and culture of doliolids

Laboratory cultures of D. gegenbauri were initiated and maintained as previously described (Walters, Gibson, et al., 2019). Briefly, cultures were maintained in filtered (~0.45 μm) seawater in 3.9L glass jars on a slowly rotating plankton wheel at ~20°C and supplemented daily with an algal mixture consisting of Isochrysis galbana, Rhodomonas sp. and the diatom Thalassiosira weissflogii. Algal concentrations were maintained at ~80 μg C L^-1. Dolioletta gegenbauri zooids were collected using a conical net with a 0.5 m opening, 202 μm mesh net, and a 4L aquarium cod end from the mid-shelf region (25–40 m isobath) of the South Atlantic Bight (SAB) between 29.6–31.2°N and 80.1–80.3°W. The net was towed from a drifting ship by slowly lowering and raising the net through the depth of the water column. The culture was initiated in May 2015 and maintained until the experiments were completed in April 2016. Over the course of this period, the culture was supplemented twice with freshly collected D. gegenbauri zooids in August and December 2015.

Doliolid feeding experiments

Figure 1 provides an overview of the design of the feeding experiments. Independent feeding experiments were conducted in January 2016 (Exp1), March 2016 (Exp2), and April 2016 (Exp3). Briefly, seawater containing natural microbial communities was collected from the same area in the SAB where doliolids were collected prior to each experiment. Prior to use, water was stored in 20L carboys in the dark at 20°C. Water used in Experiment 1 was collected 1 month prior to the experiment on December 2, 2015, water used in Experiment 2 was collected 2 days prior to the experiment on March 16, 2016, and water used in Experiment 3 was collected 8 days prior to the experiment. Water was collected from near the bottom of the water column where particle concentrations were highest. Over the period that water was collected, near bottom water temperature ranged from 17.0–22.6°C, salinity ranged from 34.2–35.8 PSU, and chlorophyll a concentration ranged from 1.04–1.14 μg L^-1. Experiments were initiated by acclimating 15–20 D. gegenbauri zooids from the culture in a 1.9L culture jar containing fresh seawater and allowed to feed for 2 h while rotating on a plankton wheel. Following the acclimation period, the zooids were transferred to a clean 1.9L jar containing fresh seawater and allowed to feed while on the plankton wheel for an additional 2 h after which time 5 zooids were sampled. Based on the gut residence times and estimated clearance rates (Gibson & Paffenhöfer, 2000), it is expected that doliolids would have cleared 250–700 mL (13%–35%) of the feeding vessel volume during the 2 h incubation. Following the feeding period, the remaining zooids were divided into 3 clean 1.9L culture jars and returned to the plankton wheel. Two of the jars contained 10 μm filtered water depleted of algal and detrital particles but would be expected to have retained natural bacterial communities. Fecal pellet samples were collected from these jars after 2 h and 24 h, respectively. To produce starved “Empty Gut” samples to characterize a core gut microbiome, an additional 5 zooids were incubated for 24 h in a third 1.9L jar containing 0.2 μm filtered seawater. An initial series of 3 pilot experiments conducted in 2012 and following the procedures described above generated the materials used to estimate bacterial concentrations in starved (EG), aged fecal pellets (FP24h), and seawater (SW) by quantitative PCR (qPCR).

Over the course of each experiment, four types of samples were collected for 16S metabarcoding characterization (Fig. 1). Seawater (SW) containing natural microbial communities was sampled at the start of the 2 h feeding period of the experiment. Prior to filtration, the water was pre-filtered through a 63 μm sieve to remove any larger aggregates and metazoans. Triplicate 500 mL samples were gently filtered onto 25 mm 0.2 μm Supor filters (PALL Life Sciences, East Hills, NY, USA). Filters were placed in sterile 2 mL cryovials and stored at −80°C until DNA was extracted. Following the 2 h feeding period, 5 zooids with full guts (FG) were sampled to characterize microbial communities associated with both the doliolid and the prey that was consumed. These samples were expected to be dominated by fresh fecal material. Following 2 h of incubation in 10 μm filtered seawater fecal pellets were collected (FP2h). This material was considered freshly released and representative of initial colonization by natural bacterial communities. Following 24 h of incubation in 10 μm filtered seawater, a second set of fecal pellet samples were collected (FP24h). This material was considered aged fecal material and representative of material with a developed mature microbial community. It was not possible to recover intact pellets after 24 h as they were highly degraded and had either broken-up completely or dissolved.

Following collection of D. gegenbauri zooid samples, zooids were stabilized prior to DNA extraction as previously described (Walters, Lamboley, et al., 2019). Briefly, zooids were anesthetized in 0.2 μm filtered seawater containing 0.4% MS‐222 (3‐aminobenzoic acid ethylester; Alfa Aesar). After rinsing each zooid three times in fresh 0.2 μm filtered seawater containing MS‐222 they were transferred to 2 ml tubes containing extraction ATL buffer from the DNeasy Blood and Tissue kit (Qiagen, Valencia CA, USA). Fecal pellets were collected by centrifugation (500 x g for 5 min) after they were rinsed 3 times in 0.2 μm filtered seawater. Pellets were transferred to 2 mL tubes containing extraction ATL buffer. All samples were stored at 4°C until DNA was extracted, usually within 24–48 h after initial collection.

DNA extraction, PCR, and Illumina sequencing

Genomic DNA from doliolid zooids was extracted and purified using the Qiagen DNeasy Blood & Tissue kit. DNA from water samples was purified using the Qiagen PowerWater kit (Qiagen, Valencia CA, USA). Manufacturer instructions were followed for both kits. Purified DNA extracts were quantified using a Qubit® 2.0 Fluorometer with the dsDNA HS assay reagents (ThermoFisher Scientific). Yields ranged from 40–254 ng (0.20–1.3 ng/μL) DNA per gonozooid and 0.26–0.49 ng/μL DNA per 100 mL of water. DNA samples were stored at −20°C until further analysis.

DNA extracts were used to amplify (in triplicate) the V4 region of the 16S rRNA gene using the primers 515F (Parada et al., 2016) and 806R (Apprill et al., 2015). Dual-index primer constructs were designed by modifying the Earth Microbiome Project (EMP) Illumina amplicon protocol as described in Schuelke et al. (2018). All primer constructs and oligo sequences have been made available on FigShare (https://doi.org/10.6084/m9.figshare.5701090). Amplification of the 16S rRNA gene was performed following the EMP protocols (Caporaso et al., 2012). PCR reactions contained 1 μL of DNA template, 0.5 μL of each primer (10 μM), 10 μL of Platinum Hot Start PCR Master Mix (2x) (Thermo Fisher), and 13 μL of molecular-grade water. Positive (ZymoBIOMICSTM Microbial Community Standard; Zymo Research, Irvine, CA) and negative (molecular-grade water, HyClone HyPure Water, GE, Healthcare Life Sciences) controls were included in all PCRs. Thermocycling profile included the following steps: 94°C for 3 min; 94°C for 45 s, 50°C for 60 s and 72°C for 90 s for 35 cycles; and 72°C for 10 min. Amplification success was evaluated with gel electrophoresis (agar 1%). Additional details on PCR conditions and purification are also provided in Schuelke et al. (2018).

Sample DNA concentrations were measured using a Qubit® 3.0 Fluorometer with the dsDNA HS (High sensitivity) Assay Kit (Thermo Fisher Scientific) and normalized prior to pooling. The DNA library was subjected to a final magnetic bead cleanup step, followed by size selection (300–700 bp range) on a BluePippin (Sage Science, Beverly, MA) to remove any nontargeted DNA. A Bioanalyzer trace was run on the size-selected pool as a quality control measure, and the library sequenced on the Illumina MiSeq Platform (2 x 300-bp paired-end run) at the UC Davis Genomics Core Facility (Schuelke et al., 2018). All wet laboratory protocols and downstream bioinformatics scripts used in this study have been deposited on GitHub (https://github.com/BikLab/doliolids).

Quantitative polymerase chain reaction (qPCR) assay Bacterial abundances associated with EG, FP24h, and SW sample types were estimated by realtime qPCR on a Bio-Rad CFX96 Real-Time PCR System (Bio-Rad Laboratories, Hercules, California) using 932F and 1062R universal 16S rRNA targeted primers (Allen et al., 2005). qPCR reactions were performed in 20 μL reactions containing a final concentration of 1X SsoFast™ EvaGreen® Supermix (Bio-Rad Laboratories, Hercules, California), 0.3 μmol of each primer, and 1 μL genomic DNA per reaction. Quantitative standard curves were generated from a 6-order of magnitude serial dilution of plasmid DNA (pDNA) containing a cloned copy of the target 16S rRNA gene (E. coli) ranging from 10¹ to 10⁷ target gene copies per reaction. qPCR cycling conditions included the following steps: 95 °C for 30 s followed by 45 cycles of denaturation (95 °C, 5 s) and annealing/extension (54.7 °C, 5 s). After cycling, product melt temperatures were evaluated from 60 to 95 °C at 0.5 °C increments for 5 s each. Samples, standards, and no-template controls were assayed in triplicate. rRNA gene copy numbers were normalized to volume. Volumes of doliolids and fecal pellets were estimated from microscopically determined size estimates assuming a barrel and spherical shape, respectively.

Bioinformatics and statistical analyses Raw Illumina data were demultiplexed using a custom script for handling dual-index barcodes and then analyzed in QIIME2 version 2020.11 (Bolyen et al., 2019). Primer and adapter sequences were trimmed (error rate of 0.1 and reads discarded when lacking adapter/primer) using the cutadapt plugin (Martin, 2011). Denoising was based on optimal parameters (forward and reverse reads truncated at 237 and 253 bp, respectively; median PHRED score of ≥30). Amplicon Sequence Variants (ASVs) were estimated using the high-resolution single-nucleotide difference DADA2 method (Callahan et al., 2016) based on default parameters and consensus chimera checking parameters. Taxonomy assignments of ASVs were obtained with the BLAST+ consensus taxonomy classifier [minimum confidence value of 0.8; (Camacho et al., 2009)] and the SILVA 138 SSURef NR99 release (Quast et al., 2013).

Our final dataset consisted of 115 samples, and except for PCR negative controls, all had high sequence depth (>2,000 reads, Table S1, Appendix S1). Preliminary analyses revealed that ASVs found in PCR controls were not shared by experimental samples (Fig. S1, Appendix S2). Still, the R package decontam (prevalence; threshold of 0.5) was used to assess the levels of contamination in the dataset (Davis et al., 2018). Only four ASVs were identified as contaminants, which were removed prior to analyses assessing patterns of microbial community variation among sample types.

Diversity estimates including Observed, Shannon H′ (Log2), Inverted Simpson (D), and Pielou’s Evenness (J’) were extracted from the filtered ASV table using the package phyloseq (McMurdie & Holmes, 2013), and compared among sample types for each experiment separately. Data normality was assessed using Shapiro–Wilk's method, and Kruskal–Wallis (K–W) tests were used to assess differences among sample types with the package FSA v0.8.24 in R version 4.1.2 (R Core Team, 2021). The Mann–Whitney U test with adjustments for p-value [BH method; (Benjamini & Hochberg, 1995)] was used for pairwise comparisons (Zar, 2010).

To visualize the similarity of microbial communities among sample types, a similarity matrix based on Bray-Curtis similarity and ASV-transformed abundances standardized by total and square root transformed values was constructed. Ordination was done by Non-parametric Multidimensional Scaling (nMDS) and Goodness-of-fit given by the stress value (Clarke, 1993). The Permutational Analysis of Variance (PERMANOVA) was used to test for significance among sample types (Anderson et al., 2008).

Differential abundance analyses were performed using the R package ALDEx2 (v1.12.0) (Fernandes et al., 2013, 2014). ASV counts were transformed using the Centered-Log Ratio method for a compositionally coherent inference and estimates. Significant differences (p < 0.05) among sample types were assessed through K–W tests at each taxonomic rank. False discovery rates (FDRs) were estimated using the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995). A heatmap (taxa: rows; samples: columns) depicting the variation among sample types of taxa differentially abundant was produced. PICRUSt2 was used to predict potential gene functions from microbial community profiles associated with sample types (Douglas et al., 2020). Differential abundance analyses on the matrices of predicted gene functions were performed and visualized as described above. Additionally, predicted gene functions were organized into distinct metabolic pathways following Yilmaz et al. (2015). In this study, all visualizations were produced with ggplot2 v.3.1.1 (Wickham, 2016) using R v.4.1.2 (R Core Team, 2021).

Raw fastq files (gz format) are uploaded in a folder. File names correspond to the sample types/treatment used in this study. See metadata for additional information regarding files.