DOI: 10.5061/dryad.kd51c5b8s
title: The microbiome of the pelagic tunicate Dolioletta gegenbauri: A potential link between the grazing and microbial food web

description: 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.
description: <strong>Collection and
 culture of doliolids</strong> Laboratory cultures of <em>D.
 gegenbauri</em> 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 <em>Isochrysis galbana</em>,
 <em>Rhodomonas</em> sp. and the diatom
 <em>Thalassiosira </em><em>weissflogii</em>. Algal
 concentrations were maintained at ~80 μg C L<sup>-1</sup>.
 <em>Dolioletta gegenbauri</em> 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
 <em>D. gegenbauri</em> zooids in August and December
 2015. <strong>Doliolid feeding
 experiments</strong> 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<sup>-1</sup>. Experiments were initiated by acclimating
 15–20 <em>D. gegenbauri</em> 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 <em>D.
 gegenbauri</em> 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. <strong>DNA extraction, PCR,
 and Illumina sequencing</strong> 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<sup>1</sup> to 10<sup>7</sup> 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 (&gt;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 &lt; 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).
description: 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.