Animals and experimental infection
All the protocols of animal handling and sampling were approved by the Animal Care and Ethics Committee, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences. All the methods were carried out in accordance with the approved protocols and relevant guidelines.
Four hundred virus-free H. diversicolor supertexta (size range between 49.73 and 58.24 mm) were bought from Xiamen, Fujian Province, and transferred to Qingdao, Shandong Province of the China country by air in April 18th, 2016. These abalones were maintained in 50 L tanks (40 abalones per tank) supplied with aerated, filtered seawater and adequate seaweed (Laminaria japonica). At the end of the two-week of acclimation period, 30 animals were selected randomly and tested negative for HaHV-1 DNA. The salinity and temperature of water were fixed at 30 ± 1 ppt and 19 ± 1 °C, respectively and half-changed daily throughout the experiment.
A viral inoculum was prepared from H. diversicolor supertexta found infected with high HaHV-1 loads, originally collected from abalone farms in the Guangdong Province in 2003. The standard protocol for the OsHV-1 inoculum preparation, described in42, was employed to prepare tissue homogenates, except that 0.22 µm-filtered natural seawater instead of artificial seawater was used in all dilution steps. Tissue homogenates for negative controls were prepared using HaHV-1-negative H. diversicolor supertexta with the same protocol. Filtered tissue homogenates were stocked at 4 °C until use. An aliquot of each tissue homogenate (200 µL) was used for HaHV-1 DNA detection and quantification by quantitative PCR (qPCR) (described below).
For the experimental infection, abalones were firstly anaesthetized with 10 g/L MgCl2, and then abalones were randomly divided into challenged and negative control groups of 180 and 70 individuals, respectively. For the challenged group, 100 μL of viral inoculum (adjusted at 104 copies of viral DNA/μL) was injected into the pedal muscle of 180 abalones: 150 of them were maintained in three 50 L tanks and used for serial sampling whereas 30 of them were placed in three 18 L tanks and used to record survival. For the negative control group, 100 μL of control homogenate was injected. A total of 40 animals were maintained in one 50 L tanks and used for serial sampling, whereas 30 were placed in three 18 L tanks and used to record survival.
Animal sampling and HaHV-1 DNA quantification
Six (2 abalones per tank) and 3 abalones were sampled at 0, 12, 24, 30, 36, 48, 60 and 72 hours post injection (hpi) from challenged and negative control groups, respectively. Four types of tissue (mantle, gill, hepatopancreas and neural tissue surrounded by some muscle) were dissected from each individual and divided in 2 pieces for DNA and RNA extraction. DNA extraction was performed from tissues and homogenates with a TIANampTM Marine Animals DNA Kit (Tiangen Biotech, China), according to the manufacturer’s protocol. The purity and quantity of the isolated DNA was determined with a Nanodrop 2000 spectrophotometer (Thermo Fischer Scientific, Germany).
HaHV-1 DNA quantification was carried out by qPCR targeting ORF66 (annotated as unknown protein) and using a protocol adapted from the World Organization for Animal Health (OIE) Manual of Diagnostic Tests for Aquatic Animals, 2017. Briefly, amplification was performed in 25 µl reactions containing 12.5 µl of 2x FastStart Essential DNA Probes Master (Roche Diagnostics, Swiss), 1 µl of each primer (10 µM), 0.5 µl of TaqMan® probes targeting the viral ORF66 (10 µM), 2 µl of template DNA and 8 µl of water. The PCR assay was performed using a Bio-Rad CFX Connect Real-Time system (Bio-Rad Laboratories, USA) and run under the following conditions: 1 cycle 95 °C for 10 min, followed by 40 cycles of amplification at 95 °C for 10 s, 60 °C for 20 s. The virus quantitation was carried out by comparison with a standard curve, which was created from a 10-fold dilution series (108−101 copies µl−1) of plasmid containing the target sequence. A qPCR negative control was carried out with 2 µl of deionized water instead of DNA sample. Each sample was tested in duplicate and it was recorded as positive if both replicates were positively amplified. We estimated the HaHV-1 infection burden of each sample as the mean genomic equivalent (GE) score (ng−1 of total DNA) of the duplicates.
RNA extraction and sequencing
Due to the excessive polysaccharide content of abalone samples, attempts to extract high quality RNA for high-throughput sequencing failed in our laboratory. Therefore, 3 mantle samples at 60 hpi were sent to Beijing Novogene Technology Co. Ltd. (China) for RNA extraction and whole transcriptome sequencing based on Illumina technology. A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations and index codes were added to attribute the reads to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads (NEB). Fragmentation was carried out using divalent cations under elevated temperature in NEB NextFirst® Strand Synthesis Reaction Buffer (5x). First strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H−) (NEB). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H (NEB). Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities (NEB). After adenylation of the 3′ ends of DNA fragments, NEBNext® Adaptor with hairpin loop structure were ligated to prepare for hybridization. To select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, USA). Then 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Subsequently, PCR was performed with Phusion High-Fidelity DNA polymerase (NEB), Universal PCR primers and Index (X) Primer (NEB). At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions and sequencing was carried out on an Illumina Hiseq platform (2 × 150 paired-end reads).
Viral genome sequencing
To reduce the noise of host DNAs, a Long-Range PCR (LR-PCR) based approach was used to enrich HaHV-1 DNA for Pacbio sequencing. Twenty one primer pairs were designed with the online version of GenoFrag based on the genome sequence of HaHV-1 Taiwan variant (GenBank accession no. KU096999)43. The length of LR-PCR products was set at 9–15 kb and overlapped the adjacent ones by 500–1500 bp. DNA was extracted from challenged samples using the Qiagen DNeasy Blood and Tissue kit (Qiagen, USA) according to the manufacturer’s protocol. PCR was performed in a 50 μl reaction volume that included 10 µL 5x PrimeSTAR GXL Buffer, 4 µL dNTP Mixture (2.5 mM each), 1 µL PrimeSTAR GXL DNA Polymerase (Takara, Japan), 1 µL each of forward and reverse primers (10 µM, listed in Supplementary Table 6), 2 µL DNA template and 31 µL of PCR-grade H2O. LR-PCR was performed using Veriti Thermal Cycler (Applied Biosystems) under the following conditions: 94 °C for 1 minutes, followed by 35 cycles of amplification (denaturation 98 °C, 10 s; annealing 50 °C, 15 s, extension 68 °C, 10 minutes) and hold at 4 °C. PCR product sizes were detected on 0.8% 1 × TAE agarose gels stained with GeneFinder™ (Zeesan Biotechnology Inc.). PCR amplicons were purified with QIAquick PCR Purification Kit (Qiagen) and quantified using Qubit Fluorometer (Life Technologies). Then, the amplicons were mixed in equimolecular proportions, and 10 µg of the mixtures were send for genome sequencing and assembly at Guangzhou Gene Denovo Biotechnology Co., Ltd. Sequencing was performed on the PacBio RS II sequencer with 10 kb SMRTBell library prepared with manufacturer’s specification (Pacific Biosciences, USA).
If not differently indicated, all the analyses were performed using CLC Genomic Workbench v.11.0 (Qiagen). Raw reads were trimmed for the presence of adaptor sequences and for quality using TrimGalore! (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), allowing a maximum of 2 ambiguous bases and a quality threshold of Q20. The 8 available Malacoherpesviridae genomes were retrieved from the NCBI database (IDs: AY509253, GQ153938, KY242785, KY271630, MG561751, KP412538, NC018874 and KU096999). High quality (HQ) reads were mapped on all Malacoherpesviridae genomes (with similarity and length mapping parameters ranging between 0.8/0.5 and 0.9/0.9, respectively) and mapped reads were labeled as ‘Malacoherpesviridae reads’ and retained for subsequent analyses. A large gap read mapping tool was employed to identify spliced reads, i.e. reads mapping on the reference genome with an internal gap.
Analysis of DNA reads
Firstly, the resulting Pacbio reads longer than 500 bp with a quality value over 0.75 were merged together into a single dataset. Then, the hierarchical genome-assembly process (HGAP) pipeline44 was used to correct for random errors in the long seed reads with a threshold of 6 Kb, by aligning shorter reads against them. Finally, de novo assembly was performed using the corrected and preassembled reads by Celera Assembler with an overlap-layout-consensus (OLC) strategy45. Since SMRT sequencing introduce very little variations of the quality throughout the reads, no quality values were used during the assembly. To validate the quality of the assembly and determine the final genome sequence, the Quiver consensus algorithm was used44.
De-novo and viral protein characterization
De-novo assembly of ‘Malacoherpesviridae reads’ was performed with the CLC assembler tool, setting word and bubble size parameters to ‘automatic’ and a minimal contig length of 500 bp; the assembled contigs were subjected to open reading frame (ORF) prediction with the NCBI ORF finder tool (https://www.ncbi.nlm.nih.gov/orffinder/), setting a minimal length of 100 codons according to the criteria described by2 and used for the annotation of all Malacoherpesviridae genomes. Briefly, overlapping ORFs and ORFs shorter than the minimal length were retained if displayed features supporting their effective existence, like conserved domain(s), a transmembrane or a signal peptide region. SignalP and HMMer46 were used to identify a signal peptide region or the presence of conserved protein domains, respectively (using the Pfam-A models, v.2947, with a cut-off E-value of 0.01).
Identification of variable positions on viral genome
Single Nucleotide Polymorphism (SNP) analysis was separately performed on the 3 mapping files (produced mapping the reads on KU096999.1 genome reference). Nucleotide changes were called ‘SNP’ if present at least in 5% of the locally aligned reads using the following parameters: minimum average quality of the five surrounding bases, 30; minimum quality of central base, 30; minimum required coverage, 50x. InDels were identified using an algorithm that first identifies positions in the mapping files with an excess of reads with unaligned ends. Once these positions and the consensus sequences of the unaligned ends were determined, the algorithm maps the consensus sequences to the reference sequence around other positions with unaligned ends. Subsequently, the consensus assembled from the DNA reads were used to validate SNPs and InDels, by mapping them on the reference genome and manually inspecting the positions of all the predicted variations.
In-silico expression analysis
To quantify the expression of viral ORFs, HQ RNA reads were mapped on the improved version of the HaHV-1 genome (KU096999.1), setting both length and similarity parameters to 0.9. Starting from the read counts per ORF, Transcripts Per Million (TPM) expression values were computed according to48.
Viral gene expression analysis
To further study the viral gene expression over the time course experiment, 37 primer pairs for qRT-PCR analysis were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Supplementary Table 7). These 37 ORFs were selected according to their expression levels of the RNA-seq data at 60 hpi. The efficiency of each primer pair was verified to be within the range of 90–110% of efficiency by constructing a standard curve from serial dilutions (Supplementary Table 7).
For qRT-PCR, total RNA was prepared in our lab using TRIzol reagents (Invitrogen, USA) in accordance with the manufacturer’s instructions. cDNA was synthesized from 2 μg of total RNA with reverse transcriptase (Takara, Japan) and random primers. qRT-PCR was performed in a total of 20 μL reaction system using TB Green™ Premix Ex Taq™ II (Takara, Japan) based on CFX Connect™ Real-Time System (Bio-Rad Laboratories, Inc.). Each reaction system contained 10 μL of TB Green Premix Ex Taq (Tli RNaseH Plus), 0.4 ul of ROX Reference Dye II, 0.8 ul of each primer at the final concentration of 400 nM each, 6 ul of distilled water, and 2 ul of cDNA dilution (1/20). qPCR reactions were performed under the following thermal cycling conditions: 1 cycle of 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 60 °C for 30 s and a melt curve step (from 65 °C, gradually increasing 0.5 °C/s to 95 °C, with acquisition data every 1 s). The relative expression levels of viral genes among samples were normalized to the abalone cytochrome c oxidase subunit I (COX I), which has been verified as reliable internal standard (unpublished data). All reactions were performed in triplicates and the data were calculated as the mean of relative mRNA expression using the 1/delta Ct method49 and were therefore reported in a scale for 0 to 1.
Phylogenetic analysis was carried out on 40 concatenated homologue ORFs retrieved from the 3 HaHV-1 genomes. Briefly, the homologue ORFs were retrieved from the Australian viral variant (NC018874), from the improved version of the Taiwanese variant herein described (KU096999.1) and from the genomic consensus obtained through PacBio sequencing. ORFs were concatenated and aligned using MUSCLE50 and a phylogenetic tree was produced based on NJ algorithm and Jukes-Cantor distance measurement, applying 1,000 bootstrap replicates.
The high-quality RNA-Seq reads and raw PacBio reads are available through the SRA database under accession numbers PRJNA47124 and PRJNA492770 respectively.