The microbiome plays a key role in animal health and is an important aspect of any natural or farmed ecosystem. Here we present the first environmental microbiome study of Ostrea edulis, as well as the first of a natural wild spawning event of any oyster species. Larval abundance was hypothesised to be correlated with specific microbial signatures. Water samples were collected throughout a natural spawning event of O. edulis at Loch Ryan, Scotland, UK. Samples were collected on 4 different dates from June to September 2019, across 8 different sampling sites on the Loch at mid, bottom, and surface levels within the water column to remove the effects of salinity and tidal fluctuations. Larval count data was obtained from these samples before full-length sequencing of the 16S rRNA gene using Oxford Nanopore Technologies. Significant microbial differences were only found between samples collected on different dates, and not at different sites or water column depths. Differences in the microbiome throughout the spawning season were driven by changes in the abundance of certain taxa, most notably those belonging to the Rhodobacteraceae family. Inverse abundance profiles of Rhodobacteraceae and Vibrio over time are also discussed. This study provides important microbial baseline data about the spawning environment of O. edulis.

Sample collection: 

Weekly water samples were collected from Loch Ryan between the 7^th June to 24^th September 2019 as described by Chapman et al, (2021). This sampling time frame was chosen to cover the native oyster spawning season as comprehensively as possible. Water samples were collected across 8 different sites in the Loch, covering 4 sites from known oyster fishing grounds off Lefnoll Point, and 4 sites from non-oyster bed habitats (see Figure 1a, Chapman et al, 2021). The location of all sample sites visited on Loch Ryan in relation to the oyster beds can be seen in Figure 1. Water was collected from the surface, mid and bottom points within the water column using the pump and trap method (see Fig.2. Chapman et al, 2021). Each sample was filtered through a 90µm mesh to capture plankton (but retain oyster larvae) within the water with 200 L of filtrate being measured from each sample. All samples were then used to conduct larval counts of both O. edulis and other bivalves such as cockles, mussels, and saddle oysters (Chapman et al, 2021). Water samples were then frozen in a -20°C freezer.

To aid in the downstream analysis, count data of total bivalves and O. edulis larvae were used to assign each individual water sample a “low”, “medium” or “high” count category for these two factors. This was done by dividing the highest count number recorded for bivalves and O. edulis by 3 in order to set threshold values for each category. The majority of count data for both bivalves and O. edulis were categorised at low, with just 4 and 2 samples belonging in the high categories for bivalve and O. edulis counts respectively.

A range of physicochemical parameters were recorded at the time of sample collection. Loggers measured water temperature, conductivity, and dissolved oxygen. The recording of daily sea surface temperature allowed for the calculation of ‘temperature sum’ as defined in Chapman et al, 2021. The tidal range and lunar cycle were also recorded for each of the sampling dates.

Molecular biology:

After the initial filtration and freezing, thawed water samples were thoroughly mixed by vortexing for 20 seconds. 15 ml of each was added to a 50 ml falcon tube containing 33 ml of absolute ethanol and 1.5 ml of 3 M sodium acetate. After 1 hour incubation at -80°C, tubes were centrifuged at 4°C, 4000 rpm for 30 minutes to pellet any cellular material and free DNA. Supernatant was poured off and ethanol was allowed to evaporate for 15 mins before 2 ml Liftons buffer with 100 µg/ml Proteinase K was added. This solution was lysed at 55°C for 1 hour, then spun down. 400 µl of supernatant was mixed with 800 µl of Qiagen AL buffer with ethanol, loaded onto a Qiagen spin column, and cleaned by following the remainder of the Qiagen blood and tissue kit tissue protocol (Qiagen, Hilden, Germany). In a 96-well plate setup, each DNA sample was then quantified using the QuantiFluor dsDNA system from Promega (Promega, Madison, US) before using as template for amplification of full-length 16S rRNA gene by PCR. This PCR was optimised by using universal 16S primers adapted from Callahan et al, (2019) to also detect bacteria of the Vibrio genus: F= AGRGTTYGATYMTGGCTCAG R= VGYTACCTTGTTACGACTT with Phusion High-Fidelity DNA polymerase from APExBIO (ApexBIO, Texas, US). Three positive controls were included by the use of Streptococcus bacterial DNA, alongside four negative controls of nuclease free water. The PCR reaction was run with an initial denaturation step at 98°C for 30 seconds, followed by 34 cycles of 98°C for 10 seconds, 65°C for 30 seconds, 72°C for 50 seconds in a PTC-225 Peltier MJ Research Thermocycler (Marshall Scientific, Hampton, US). The subsequent PCR product was then visualised via gel electrophoresis to check for DNA fragments of 1466 base pairs before purifying with the High Prep PCR clean-up system (MagBio, Maryland, US). Each DNA sample was then re-quantified using Qubit dsDNA BR Assay kit (Invitrogen, Massachusetts, US) and standardised to ~200 fmol through dilution with nuclease-free water. The SQK-LSK110 sequencing ligation kit and EXP-NBD196 Native Barcode Expansion kit (Oxford Nanopore Technologies, Oxford, UK), were then used to perform library preparation with metabarcoding on each of the samples. The pooled library was then loaded across 4 separate ONT MinION R10 flow cells for 48 hours, obtaining a total of 98.681 GB of data.

Data handling and analysis:

High accuracy basecalling was performed using Guppy v.6.1.7 (ONT) and reads were quality controlled by inspecting the average quality and length with Nanoplot version 1.4.0 (De coster et al, 2018). Trimming of primers from reads was then performed using Cutadapt v.3.4 (Martin, 2011), with a maximum error rate of 0.15 and a minimum overlap of 20 base pairs (bp). Reads shorter than 1,300 bp and longer than 1,600 bp were also discarded. Porechop v0.2.4 (Wick et al, 2017) was then used to trim reads of chimeras and adapter sequences. Microbial taxonomy was then identified using Kraken2 v2.1.1 (Wood and Salzberg, 2014), where 2,112,739 reads were aligned against the SILVA138 database (Wood, Lu, and Langmead, 2019).

Downstream data analysis was carried out in R studio v4.2.0 (RStudio Team, 2022) using the Phyloseq (McMurdie and Holmes, 2013) and ggplot2 (Wickham, 2014) packages. This included the calculation of observed and Choa1 alpha diversity measures which were then tested for significant differences with a Kruskal-Wallis test. Control and outlying samples were identified from non-metric multidimensional scaling (NMDS) analysis and excluded, leaving 85 samples in the results presented. Use of the Vegan package (Oksanen, 2012) in R Studio also allowed for a PERMANOVA statistical test to be performed, and weighted phylogenetic tress (see Figures 3 and 6), were created using the Pavian web application tool (Breitwieser and Salzberg, 2020).