2.1 Pacific oyster family production and maintenance
Ten biparental Pacific oyster (Crassostrea gigas) families were produced in March 2019 at the Cawthron Institute’s hatchery in Nelson, New Zealand (41°11'33.3"S, 173°21'37.8"E). Broodstock selected for spawning originated from the Mahurangi harbour (Warkworth, New-Zealand; 36°25'28.16"S, 174°41'32.36"E). Parents used to produce resilient families have been selected for three generations and were survivors from full-sib families that had been exposed to an on-farm virus challenge and selected based on their high survival rates in the field. Parents used to produce susceptible families were derived from i) a subset of the families that showed poor survival during the on-farm challenge; and ii) specific males that were reared on an uninfected farm on the South Island, and were, therefore, expected to be naïve and highly susceptible to the virus [12]. Cryopreserved sperm from these individual naïve males were thawed according to Adams et al. [33].
Twenty million eggs from each female were fertilized with either fresh sperm or thawed sperm from a single male. Resulting embryos from a single cross were then incubated separately for 24h according to [34]. Individual families were reared using a high-density larval rearing flow-through system in the hatchery for 5 weeks according to a modified protocol from Ragg et al. [35] .
Oyster families were then further on-grown for 5 months in a common flow-through water tank but kept separate in family-specific upwelling tanks. All families were reared under common conditions: seawater temperature ranged from 9°C to 18°C, and algal food was continuously supplied to provide an optimal growing environment for the spat. Husbandry treatments such as grading and biomass adjustments, were conducted simultaneously between families. Despite these measures, mean liveweight of spat at the beginning of the trial varied among families (Supplementary Table 1). In the nursery tank, oyster families were continuously supplied with UV-sterilized seawater (80 mJ cm-2) and maintained under strict biosecurity management to ensure that no pathogen would interfere with later experiments. The “pathogen free” status of the experimental oysters was confirmed prior to the initiation of the experiment: indeed, no OsHV-1 DNA was detected using qPCR (n=100; [36]). Finally, no significant spat mortality was observed prior to the start of the experiment.
2.2. Experimental design
2.2.1. Acclimation of oysters
On the 4th September 2019, virus-free oysters were randomly collected from their family-specific upwelling tanks and transferred to the experimental challenge facility at the Cawthron Institute (Nelson, New-Zealand; 41°16'16.7"S, 173°17'36.3"E) for acclimation.
Experimental infection protocols consisted of a water transfer between tanks containing C. gigas oysters carrying the disease (referred as ‘donors’) and tanks containing pathogen-free oysters (referred as ‘recipients’) adapted from [16,18,37]. Two pools of 5,000 spat each (6-month-old oysters, mean individual weight of ~ 1.2 g), consisting of a mixture of oysters from our susceptible families, were used as “donors”, and placed in two 300L tanks; one tank for the control donors (which were to be injected with artificial seawater) and one tank for the pathogen donors (which were to be injected with OsHV-1 suspension) (Fig. 1). Tanks were continuously supplied with 1 µm-filtered and UV-sterilized seawater (FSW), preheated through a Digiheat inline heater (Waterco, Auckland, New Zealand), and flow rates were maintained at 1.5L min-1 per tank. Donors were maintained in these conditions for 2 weeks until their injection.
In parallel, family-specific mesh bags containing 195 oysters (6-month-old spat, mean individual weight of ~ 1.2 g) were prepared for each of the 10 families and used as “recipients”. One bag per family was randomly suspended (i.e., 10 families x 1 bag = 10 bags per tank) into six 100L tanks. Three replicate recipient tanks were set up for the pathogen group and three replicates for the control group (Fig. 1 A, B). Total biomass per recipient tank (corresponding to 1,950 spat) was 2.4 kg ± 0.2. Flow rates were adjusted to 0.5L min-1 for each tank by means of a valve so combined flow rate for the triplicate tanks was 1.5L min-1. Recipient oysters were acclimated under these conditions for 15 days until viral inoculation. No mortality was observed during this time.
To thoroughly maintain seawater temperature at 20.8°C (± 0.8), oxygen saturation above 90% and seawater well homogenized, all tanks were equipped with an aquarium immersion heater (Eheim thermocontrol 200W), light aeration, and a circulation pump (Hailea, Low water level pump DS-700). All oysters (donors and recipients) were fed ad libitum with a bispecific mixture of hatchery-grown Chaetoceros muelleri (CS-176) and Tisochrysis lutea (CS-177) continuously supplied using a peristaltic pump. Algal background concentration was maintained between 3 to 10 µg L-1 Chlα (equivalent to 5 to 20 cells of Tiso per µL).
2.2.2 Viral suspension
The OsHV-1 suspension stock was produced in April 2014 as described in Camara et al. [12]. Briefly, tissue from high virus load oysters was homogenized, cell debris was removed by centrifugation, and the supernatant was purified by serial filtrations to 0.22 μm. Finally, a cryoprotectant solution (10% glycerol and 10% fetal calf serum final concentration) was added and the resulting suspension slowly frozen and stored at −80 °C.
2.2.3 Infection and sampling procedures
On the 17th September 2019, pathogen-donor oysters were myorelaxed in hexahydrate MgCl2 (30 g L-1) according to [38] until valve opening. Concurrently, cryopreserved virus stock suspension was thawed in a 22°C water bath for 10 minutes, and diluted 1/5 in sterile artificial seawater (SSW). Pathogen donors were then injected in the adductor muscle using a 26-gauge needle attached to a multi-dispensing hand pipette, with 20 μl of viral suspension (1.76 x 105 copies of OsHV-1 per injection), while “control donors” were injected with the same volume of sterile SSW. Oyster donors (pathogen and control) were then held at 22°C (± 0.5) for 70 h in their respective tanks, in static conditions (i.e., no water renewal) to produce infectious or control water. During the incubation period, survival of donor oysters was monitored daily, and dead animals were immediately removed from the tanks.
On the 20th September, infectious or control water was transferred to the respective recipient tanks (3 tanks OsHV-1 challenged, 3 tanks control). Feeding was stopped and the recipient tanks were left in static conditions at 22°C (± 0.5). After 48 hours, new FSW was gradually reintroduced to the recipient tanks at a flow rate of 0.5 L min-1. Recipients were then inspected daily and dead animals, characterized by failure to close their valves, were immediately removed. Survival of recipients was monitored for 14 days or 336 hours post infection (hpi). To avoid accidental releases of OsHV-1 in the environment, continuous chlorination (200 ppm for 1h) combined with UV sterilization (80 mJ cm-2) were applied to the effluent water.
Five live recipient oysters were randomly sampled from each family-specific bag at 0, 2, 6, 12, 24 hpi, and every 24h until 120 hpi. A final sampling of 5 live recipient oysters was conducted at the end of the challenge, at 336 hpi (Fig. 1). Whole tissues were removed from the shells, dried by dabbing on a paper tissue, flash frozen in liquid nitrogen and later reduced to powder (Mixer Mill MM400, Retsch GmbH, Germany) and stored at -80°C for OsHV-1 DNA, viral gene expression and metabarcoding analyses. In addition, three live recipients were randomly collected from each family-specific bag at 0 and 72 hpi. Their tissues were carefully dissected, placed in histological cassettes, fixed in 4% formalin for 48h and stored in 70% ethanol for later histopathological analysis.
One liter of seawater was collected from each donor tank (n = 2) at 0, 48 and 72 hpi, while 1L of seawater from each recipient tank (n = 6) was sampled at 0, 24, 48, 72, 96, 120, 168, 240, and 288 hpi (Fig. 1). Water was collected using sterile (autoclaved) glass bottles, and immediately filtered through a sterile 47 mm cellulose membrane filter of 0.22 µm pore size to isolate bacteria and viruses. Membrane filters were flash frozen and stored at -80°C until DNA extraction.
2.3 Biometric analyses and water quality measurements
Each family-specific mesh bag was weighed on the 4th of September prior to acclimation (ST 1), on T0 (prior to viral infection), and at the end of the challenge on the 4th of October (T336h).
Water quality parameters were measured throughout the challenge by means of a YSI ProSolo digital meter (Xylem Inc., Yellow Springs, OH, USA) for temperature and dissolved oxygen (D.O), a handheld Testo 206 pH meter for pH, and a FluoroSense handheld fluorometer (Turner Designs, San Jose, CA, USA) for microalgal background levels.
2.4 DNA extraction and OsHV-1 DNA quantification
Total DNA was extracted from powdered donors (0 hpi) and recipient oysters (0 and 120 hpi), and from seawater samples from donor tanks (0, 48 and 72 hpi) and recipient tanks (0, 24, 48, 72, 96, 120, 168, 240, 288, 336 hpi) using Blood and tissues kit (QIAGEN) according to the manufacturer’s protocol. Four blank DNA extractions were included in order to test for potential bacterial contamination of the DNA extraction kit and/or reagents.
Level of OsHV-1 DNA was quantified using real time PCR in water from both donors (0, 48 and 72 hpi) and recipient tanks (0, 24, 48, 72, 96, 120, 168, 240 hpi) (Fig. 1). Real time PCR was carried out in 20 µl reaction mixture consisting of 10 µl of SsoFast Tm EvaGreen Supermix (Biorad), 1 µl of each primer (OsHV-1 BF and OsHV-1 B4, 10 µM), 6 µl of water and 2 µl DNA sample. PCR amplification was performed using Rotor-Gene R Q (Qiagen) following: 1 cycle pre-incubation at 98°C for 2 min, 40 cycles of amplification at 98°C for 15 s and 58°C for 20 s; melting temperature curve Ramping from 72°C to 95°C, rising by 1 degree each 5 second. Samples were analysed in triplicate, and 3 controls were carried out: a negative control which contained PCR reaction mixture without the target, an extraction control and a positive control which holds DNA target(s). The standard curve was prepared using serial dilutions of chromosomal DNA from Sydney University (Sydney School of Veterinary Science, Infectious Disease Laboratory, Farm Animal Health, Australia), to determine OsHV-1 DNA concentration.
2.5 Microbiota sequencing
Polymerase chain reactions (PCR) were performed on some of the samples extracted as described earlier (n = 120 from oyster tissue.; n = 54 from water, plus blank controls). Bacterial communities were amplified using the 16S rRNA gene (v3-v4 region) with the primer set 341F: 5’-CCT ACG GGN GGC WGC AG-3’[39] and 805R: 5’-GAC TAC HVG GGT ATC TAA TCC-3’ [40]. These primers were modified to include IlluminaTM overhang adaptors (forward: 5'-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG-3' and reverse: 5'- GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA G-3'), as described in Kozich et al. [41].
Polymerase chain reactions were carried out in 50 μL reaction volumes containing 25 μL MyFiTM 2× PCR supermix (Bioline, London, UK), 19 μL of nuclease-free H20, 0.20 μM of modified Illumina overhang adaptor primers, 2.0 μM of both blocking primers, and 1 μL of template DNA. Thermocycling conditions were 95°C for 2 min, followed by 39 cycles of 94°C for 20 s, 52°C for 20 s, 72°C for 30s, with a final extension step at 72°C for 5 min. Amplicons were purified using the SequalPrepTM normalization plate kit (Applied Biosystems, CA, USA) resulting in an equimolar concentration of ~1 ng μL-1, all according to the manufacturer’s instructions. Purified amplicons were individually indexed using the NexteraTM DNA library Prep Kit (Illumina, California, USA) and paired end sequenced on a MiSeqTM Illumina with the 2x250 bp v2 kit at Auckland Genomics (Auckland, New-Zealand). Raw sequences are publicly available from the NCBI database under project number PRJNA832870.
2.6 Viral gene expression
Viral gene expression was quantified in recipient oysters from the ten families sampled at 6, 12, 24, 48, 96 and 120 hpi. RNA was extracted from 30 mg of powdered oysters using the Direct-Zol RNA Miniprep kit (Zymo research) according to the manufacturer’s protocol. Samples were then treated with DNAse I (TURBOTM DNase, Invitrogen) to remove genomic DNA. To confirm the absence of DNA in the sampled RNA, a 16S PCR assay was performed on each RNA sample after DNase treatment and gels were run. The quality and purity of the isolated RNA in all samples were checked using a NanoPhotometer (Implen, Munich, Germany). DNAse-treated RNA was transcribed into cDNA, using the Super-Script III reverse transcriptase (Life Technologies, CA, USA). Droplet digital PCR (ddPCR) was conducted in an automated droplet generator (QX200 Droplet Digital PCR SystemTM, BioRad) to determine the expression of three viral genes ORF 27, ORF 38 and ORF 87 selected from among the 39 ORFs described by Segarra et al. [42]. These ORFs encode for different protein functions and expressed differently during an OsHV-1 replication cycle [42]. Each ddPCR reaction included 1 µl of 3 µM of each primer, 10 µl ddPCR Supermix for Evagreen (BioRad), 1 µl DNA and 7 µl sterile water for a total reaction volume of 21 µl. The BioRad QX200 droplet generator partitioned each reaction mixture into nano droplets by combining 20 μl of the reaction mixture with 70 μl of BioRad droplet oil. After processing, this resulted in a total nanodroplet volume of 40 μl, which was transferred to a PCR plate for amplification using the following cycling protocol: hold at 95°C for 5 s, 45 cycles of 95°C for 30 s, 60°C 1 s, and a final enzyme deactivation step at 98°C for 10 min. The plate was then analysed on the QX200 instrument. For each ddPCR plate run, at least one negative control (RNA/DNA-free water; Life Technologies), and one positive control (Gblock for each ORF tested, diluted 1/10,000) were included.
2.7 Histological analyses
Recipient oysters were collected for histopathological analysis at 0 and 72 hpi. Three individuals per family (10 families at 0 and exposed oysters at 72hpi and 3 families for non-exposed oysters at 72hpi) were sampled per tank (3 tanks). Following storage in 70% ethanol, whole oysters were embedded in paraffin wax, before being serially sectioned to a thickness of 5 μm using a microtome. Tissue sections were collected on polylysine-coated slides and stained with Giemsa (performed by Gribbles Veterinary pathology, New Zealand), which contains a mixture of azure and eosin that variably stain the basic components of the cell pink/purple (e.g cytoplasm, granules) and methylene blue, which stains the acidic components of the cells blue (e.g nucleus). The presence of histopathological features was assessed in the mantle and digestive gland of oysters using a light microscope and up to 1000X magnification (Olympus BX53 microscope with a DP22 digital camera). Pathological features were categorised as ‘0’ when absent and ‘1’ when present and were converted to a percent presence of a given feature per family, keeping tank replication (n=3).
2.8 Bioinformatic analysis
Sequence data was demultiplexed using the MiSeq Reporter (v.2) and primers removed using CUTADAPT (Version 2.6; [43]) allowing for no indels and a minimum overlap of 17 base-pairs (bp). Sequences were quality filtered (maxN=0, maxEE=c(2,2), trunQ=2), denoised, paired end merged (minOverlap=10) and chimera filtered (method=consensus) using the ‘DADA2’ R program (version 1.14; [44]). Prior to quality filtering, forward and reverse reads were truncated at 226 and 220 bp on the 5’ end, respectively, to remove the lower quality section. Amplicon sequence variants (ASVs) were taxonomically identified using the RDP Naïve Bayesian Classifier algorithm [45] (implemented in DADA2 using the SILVA ribosomal RNA gene database ([46]; version 138; https://benjjneb.github.io/dada2/training.html). Unassigned ASVs and those identified as non-bacterial were discarded. Additionally, potential contaminant reads were identified and removed with the MicroDecon R package ([47]; version 1.0.2) rare ASVs (ASVs for which the sum of read did not exceed more than 2 reads in at least 3 samples) were discarded. Sequencing depth per sample was visualized with the “rarecurve” function of the “vegan” R package (version 2.5.7; [48]), and samples with less than 10,000 reads discarded to ensure that samples used in downstream analyses had sufficient sequencing depth to recover most of the diversity.
2.9 Statistical analyses
Survival functions were computed according to Kaplan and Meier (1958) using RStudio 4.1.0 and the “Survival” R package (V 3.2-13, [49]). Survival time was measured in hours from the injection for (i) donors (Supplementary Data 1) or (ii) from the onset of infection for the recipients (Fig. 2). The data were read as the number of dead oysters within each bag per tank at each count. Survival time curves were compared using the cox regression model [50] after adjusting for (i) injection (OsHV-1 or SSW) for donors, or (ii) for family (F1 to F10), and infection level (OsHV-1 or control) for recipient oysters, considering the random effect of the tanks and bags. The survival of control recipient oysters was not included in the statistical model because it was 100%. The proportionality of hazards (PH) was checked based on Schoenfeld residuals [51].
Mixed-design, time-repeated ANOVAs were performed to assess differences in (i) OsHV-1 DNA load in water of recipient tanks according to Family (ten levels) and time (eight levels) and ii) percent presence of a pathological feature depending on family, time, and exposure to OsHV-1. The replication unit was the tank in which the ten families were maintained. All mutual interactions among factors were tested, and Tukey’s honestly significant difference test was used as a post-hoc test. The normality of residuals and homogeneity of variances were graphically checked. Statistical analyses were performed in R studio, version 4.1.0 (R; https://www.R-project.org/). For viral gene expression, heatmaps were constructed using Multiple Experiment Viewer software ([52] http://mev.tm4.org/#/datasets/upload).
Microbial taxonomic composition in seawater and in oysters was investigated and visualised at phylum and family levels using bar plots and the ‘ggplot2’ R package (version 3.3.5; [53]). Alpha-diversity metrics such as ASV richness, Shannon and Simpson indexes were computed with the ‘Phyloseq’ R package (version 1.34.0; [54]) and visualised with line plots using ’ggplot2’. The effect of treatment, family type (vulnerable vs resistant vs highly resistant), collection date and their interactions on alpha-diversity metrics of oyster microbiota was investigated with linear mixed-effects regressions (LMER) using the ‘lme4’ R package (version 1.1.27; [55]). Similarly, the effect of treatment, time after start of experiment, and their interactions, on alpha-diversity metrics of seawater microbial communities was investigated with a LMER.
The effect of treatment on the oyster microbial community composition and structure was visualised with a principal component analysis (PCA) and tested with a permutational analysis of variance (PERMANOVA) using the ‘vegan’ R package (version 2.5.7; [48]) with the following parameters: adonis2 (Bray-Curtis distance matrix of read abundance table transformed to relative abundance ~ infection* family type + Time * family type, blocks = tank, permutations = 999, method = "bray", by = "terms"). Differences in core oyster microbiome between vulnerable and highly resistant families were investigated with heatmaps using the plot core function of the ‘microbiome’ R package (version 1.13.8; [56]) and using the ‘ComplexHeatmap’ R package (version 2.6.2;[57]).
Oyster bacterial genera associated with oyster mortality rate, ORF and viral load were investigated with Pearson correlation based on centered-log ratio transformed read abundance to account for the compositional nature of the data and visualised with a heatmap using the ‘ggcorrplot’ R package (version 0.1.3; [58]). In addition, interactions between these bacteria were assessed with the CoNet Cytoscape plugin (version 1.1.1; [59]) using Pearson correlations and the Bonferroni multiple-test correction. Variables that could best predict mortality rate, including oyster bacteria ASVs transformed to centered-log ratio, viral load and ORF concentration were identified using a Random Forest model trained (parameters: trControl = “repeatedcv”, number = 5, repeats = 3, search=”grid”; tunegrid = expand.grid(.try=c(1:8))) using the ‘caret’ R package (version 6.0.88; [60]) and visualised with bar plots using ggplot2.