ORIGINAL RESEARCH published: 25 May 2020 doi: 10.3389/fmicb.2020.00874 Whole Genome Sequencing of Hepatitis A Virus Using a PCR-Free Single-Molecule Nanopore Sequencing Approach Frederico M. Batista 1* , Tina Stapleton 1 , James A. Lowther 1 , Vera G. Fonseca 1 , Rebecca Shaw 1 , Christopher Pond 1 , David I. Walker 1 , Ronny van Aerle 1 and Jaime Martinez-Urtaza 1,2 1 International Centre of Excellence for Aquatic Animal Health, Centre for Environment Fisheries and Aquaculture Science (CEFAS), Weymouth, Dorset, United Kingdom, 2 Department of Genetics and Microbiology, Facultat de Biociències – Edifici C, Campus Universitat Autònoma de Barcelona (UAB), Barcelona, Spain Hepatitis A virus (HAV) is one of the most common causes of acute viral hepatitis in Edited by: John R. Battista, humans. Although HAV has a relatively small genome, there are several factors limiting Louisiana State University, whole genome sequencing such as PCR amplification artefacts and ambiguities in de United States novo assembly. The recently developed Oxford Nanopore technologies (ONT) allows Reviewed by: Gloria Sánchez Moragas, single-molecule sequencing of long-size fragments of DNA or RNA using PCR-free Institute of Agrochemistry and Food strategies. We have sequenced the whole genome of HAV using a PCR-free approach Technology (IATA), Spain by direct reverse-transcribed sequencing. We were able to sequence HAV cDNA and Yutaka Kawarabayasi, National Institute of Advanced obtain reads over 7 kilobases in length containing almost the whole genome of the Industrial Science and Technology virus. The comparison of these raw long nanopore reads with the HAV reference wild (AIST), Japan type revealed a nucleotide sequence identity between 81.1 and 96.6%. By de novo *Correspondence: Frederico M. Batista assembly of all HAV reads we obtained a consensus sequence of 7362 bases, with a frederico.batista@cefas.co.uk nucleotide sequence identity of 99.0% with the genome of the HAV strain pHM175/18f. When the assembly was performed using as reference the HAV strain pHM175/18f a Specialty section: This article was submitted to consensus with a sequence similarity of 99.8 % was obtained. We have also used an Evolutionary and Genomic ONT amplicon-based assay to sequence two fragments of the VP3 and VP1 regions Microbiology, a section of the journal which showed a sequence similarity of 100% with matching regions of the consensus Frontiers in Microbiology sequence obtained using the direct cDNA sequencing approach. This study showed Received: 03 November 2019 the applicability of ONT sequencing technologies to obtain the whole genome of HAV Accepted: 14 April 2020 by direct cDNA nanopore sequencing, highlighting the utility of this PCR-free approach Published: 25 May 2020 for HAV characterization and potentially other viruses of the Picornaviridae family. Citation: Batista FM, Stapleton T, Keywords: hepatitis A virus, Picornaviridae, genome, nanopore sequencing, metagenomics, SNV Lowther JA, Fonseca VG, Shaw R, Pond C, Walker DI, van Aerle R and Martinez-Urtaza J (2020) Whole Genome Sequencing of Hepatitis INTRODUCTION A Virus Using a PCR-Free Single-Molecule Nanopore Hepatitis A virus (HAV) is one of the most common causes of acute viral hepatitis in Sequencing Approach. humans. According to the World Health Organization (WHO, 2017), it is estimated that Front. Microbiol. 11:874. HAV caused approximately 11,000 deaths worldwide in 2015. This virus is commonly doi: 10.3389/fmicb.2020.00874 transmitted through the faecal-oral route via exposure to contaminated food and water Frontiers in Microbiology | www.frontiersin.org 1 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing and thus considered one of the most severe food-borne information on HAV genetic diversity in a very timely fashion. viruses (Sattar et al., 2000). HAV, which belongs to the family The main objective of this study was to develop a PCR-free Picornaviridae, is a single-stranded RNA virus with a genome of approach using long-reads to sequence the whole genome of the 7.5 kilobases of positive-strand polarity with a poly(A) stretch hepatitis A virus. at the 30 end (Daijogo and Semler, 2011). The virus replicates primarily in liver cells. The development of suitable, whole-genome sequencing MATERIALS AND METHODS (WGS) methods for foodborne viruses such as HAV is key to gain a better understanding about transmission events, Sample Preparation pathogenesis and disease-related effects (Quinones-Mateu et al., The cytopathogenic pHM175 strain of HAV used in this study 2014). Although foodborne viruses have relatively small genomes, was originally provided in 2005 by the group of Albert Bosch at there are several factors limiting the sequencing of the entire the University of Barcelona, Spain, and has since been maintained genome. Some of the limitations to obtain complete genomes through several passages at our laboratory. HAV was cultured at using high-throughput sequencing (HTS) are mainly due 37◦ C in a 5% CO2 atmosphere in fetal rhesus kidney (FRhK-4) to the shorter reads produced. Short reads not only can cell monolayers in 25 cm2 flasks using M2279 Minimum Essential impair the detection of different viral haplotypes (Takeda Medium Eagle with Earle’s salts and sodium bicarbonate (1% et al., 2019) and introduce ambiguities in de novo assemblies penicillin & streptomycin, 1% 200 mM L-glutamine solution, 1% (Schatz et al., 2010) but are also unable to resolve complex MEM non-essential amino acids (100 times), 10% for growth or regions of the genome (Somerville et al., 2019; Treangen 2% for maintenance of fetal bovine serum). A non-inoculated and Salzberg, 2011). Additionally, the PCR amplification steps flask was used as a control. At day 12, when considerable involved during library preparation can introduce considerable cytopathic effects were visible, flasks were frozen at −30◦ C. The bias, for example selective amplification and chimeras (Ardui viral particles were then harvested by thawing the frozen flask et al., 2018; Ng et al., 2018). Overall HTS can detect variant and tapping firmly against the heel of the hand as soon as frequency but it cannot accurately map haplotype differences partially defrosted, to detach the cells from the flask surface. and consequently viral clonal populations will never be assessed The content of the flasks was transferred to 50 mL centrifuge completely, namely the resistant and less abundant ones tubes and centrifuged for 10 min at 3000 × g to pellet cell (Posada-Cespedes et al., 2017). debris. The supernatant was split into 1 mL aliquots, and retained With the advent of ‘Third generation’ sequencing approaches at −80◦ C. such as the recently developed nanopore sequencing methods by Oxford Nanopore Technologies (ONT), it is now possible to HAV Quantification sequence single-molecules of particularly long-size fragments A one-step quantitative reverse transcription PCR (RT-qPCR) of DNA or RNA using amplicon-based or PCR-free strategies method was used to quantify HAV where RNA from frozen (Keller et al., 2018). Nanopore sequencing provides genome- FRhK cells infected with HAV was extracted as described by length reads in real time, generating an integral vison of the Tan et al. (2019). Briefly, 250 µL aliquots of viral supernatant genome that includes all possible mutations within a single virus were pipetted into 10 individual tubes containing 750 µL of particle (McNaughton et al., 2019). The particularly long reads TRIzol LS reagent (Invitrogen) each. After 10 min incubation generated by nanopore sequencing allow the complete genome at room temperature, 200 µL of chloroform was added and to be sequenced from single viral clones and consequently the mixture was shaken vigorously by hand for 30 s. This make possible the detection of undescribed or rare variants was followed by a 10 min incubation at room temperature (Takeda et al., 2019). The use of nanopore sequence data and centrifugation at 12,000 × g for 15 min. The aqueous will increase, without precedent our ability to characterize phase (colorless upper phase) was transferred into a 1.5 mL mutant haplotypes, which reflect the immense genetic diversity microcentrifuge tube, where equal volumes (ca. 500 µL) of of viral related sequences (aka quasispecies) within infected isopropanol and 2 µL of glycoblue precipitant were added. cells and host organisms (Andino and Domingo, 2015). The mixture was incubated overnight at −80◦ C. Samples were This will be important to further extend our understanding centrifuged at 12,000 × g for 15 min and the supernatant of sequence variants underlying phenotype and disease. was discarded. RNA was washed twice with 70% ethanol and Recently nanopore sequencing was used for whole-genome resuspended in 30 µL of nuclease-free water. For the one-Step sequencing of the Hepatitis C virus (Takeda et al., 2019; RT-qPCR, we used the primer pair HAV68 (50 -TCA CCG CCG Ueda, 2019) and Hepatitis B virus (McNaughton et al., TTT GCC TAG-30 ) and HAV240 (50 -GGA GAG CCC TGG AAG 2019). Nanopore sequencing was also used to serotype the AAA G-30 ) and the HAV150 MGB probe (50 -CCT GAA CCT Foot-and-mouth disease virus, which also belongs to the GCA GGA ATT AA-30 ) (Costafreda et al., 2006). One-Step RT- Picornaviridae family (Hansen et al., 2019). A major element qPCR PCR reactions were performed in 25 µL final volume of disease control and prevention strategies by connecting composed of 5 µL of RNA Ultrasense 5X Reaction Mix (Life human, animal and environmental health is recognized by Technologies), 1.25 µL of RNA Ultrasense enzyme mix, 0.05 µL the One Health (OH) concept (WHO, 2017). Combining of ROX reference dye, 0.25 µL of HAV68 primer (50 µM), demands of OH with third generation sequencing, such as 0.45 µL of HAV240 primer (50 µM), 1.25 µL of HAV150 probe nanopore sequencing will provide unprecedent in-depth (50 µM), 11.75 µL of distilled water, and 5.0 µL of template Frontiers in Microbiology | www.frontiersin.org 2 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing RNA. Amplification was performed by an initial incubation at second library 100 ng of oyster tissue RNA was used. For library 55◦ C for 60 min, followed by one cycle at 95◦ C for 5 min preparation, the first step started with cDNA synthesis using and 45 cycles of 15 s at 95◦ C, 1 min at 60◦ C, and 1 min at reverse transcription and strand-switching. Briefly, 100 ng of 65◦ C. Amplification was done on a QuantStudio 3 real-time PCR RNA was resuspended in 7.5 µl of RNase-free water and mixed machine (Applied Biosystems). For quantification purposes, a with 2.5 µl of the VN primer (ONT; 2 µM) that targets the standard curve was prepared by using serial dilutions of linear poly-A tail and 1 µl of dNTPs (10 mM total concentration). The dsDNA carrying the HAV target sequence (ranging between solution was mixed by gently flicking the tube, incubated for 1x105 to 1x101 copies/µL). 5 min at 65◦ C and immediately snap cooled on a freezer block. In a separate tube the following reaction mix was prepared: RNA Extraction for Sequencing 4 µl of 5X RT buffer (ONT), 1 µl of RNaseOUT (40 U/µl, Life For the preparation of the sequencing libraries, RNA was technologies), 2 µl of strand-switching Primer (10 µM, ONT) extracted from FRhK cells infected with HAV using TRizol LS and and 1 µl of RNase-free water. The mixture was mixed by gently the Zymo RNA clean & concentrator-5 (Zymo Research). Due flicking the tube then added to the snap-cooled RNA prepared to RNA limited yields we have also extracted RNA from Pacific above and incubated for 2 min at 42◦ C. One µl of Maxima oyster Crassostrea gigas gill and digestive gland tissues using the H Minus Reverse transcriptase (Life Technologies) was added same protocol to have sufficient RNA sample to sequence using and the mixture was incubated for 90 min at 42◦ C followed ONT. To prevent RNA shearing, samples were not vortexed by heat inactivation at 85◦ C for 5 min and finally held at 4◦ C. throughout the whole RNA extraction procedure. Briefly, 10 RNA was then degraded by adding 1 µl of RNase Cocktail separate aliquots of 250 µL of viral supernatant were added Enzyme Mix (Thermo Fisher) followed by incubation at 37◦ C to 750 µL of TRIzol LS reagent (Invitrogen) each (a total of for 10 min. The cDNA was then purified using AMPure XP 2500 µL of viral supernatant was processed). Oyster RNA was beads as described in the ONT protocol and eluted in 20 µl chosen since oysters are a common food matrix analyzed in of nuclease-free water. The second strand was synthetized in a our lab for the detection of foodborne viruses. For oyster tissue, 50 µl reaction composed of 25 µl of 2x LongAmp Taq Master approximately 30 mg of gill tissue was added to 1000 µL of mix (New England Biolabs), 2 µl of PR2 primer (ONT), 3 µl TRIzol LS reagent and homogenized using a disposable pestle. of nuclease free-water and 20 µl of reverse-transcribed sample The mixtures were then incubated at room temperature for from above. The mixture was incubated at 94◦ C for 1 min, 50◦ C 10 min. In a fume hood, 200 µL of chloroform were added for 1 min, 65◦ C for 15 min and hold at 4◦ C. The cDNA was to each sample and samples were vigorously shaken by hand then purified using AMPure XP beads as described in the ONT for 30 s. Following 10 min incubation at room temperature, protocol and eluted in 21 µl of nuclease-free water. In the final samples were centrifuged at 12,000 × g for 15 min and the step of the library preparation, the ends of the cDNA fragment aqueous phase was transferred into 1.5 mL tubes. One volume were repaired to create blunt ends and dA-tails were added. This of absolute ethanol was added and samples were carefully was performed by mixing the following reagents in a 0.2 ml PCR mixed by pipetting. The mixture was then processed using tube: 20 µl of the cDNA prepared above, 30 µl of nuclease-free the Zymo RNA clean & concentrator-5 Kit as described by water, 7 µl of Ultra II End-prep reaction buffer (New England the manufacturer. At the final step, the RNA was eluted in Biolabs) and 3 µl of Ultra II End-prep mix (New England 10 µL of DNase/RNase-free water. The quantity and quality Biolabs). After gentle mixing by pipetting, the reaction mix was of the RNA extracted from FRhK cells infected with HAV and incubated in a thermal cycler for 5 min at 20◦ C then 5 min at oyster gill and digestive gland tissues were assessed using the 65◦ C. The cDNA was then purified using AMPure XP beads as Qubit RNA HS Assay (Invitrogen) and the Nanodrop (Thermo described in the ONT protocol and resuspended in 22.5 µl of Fisher Scientific). The absence of HAV RNA in oyster gill nuclease-free water. Sample barcoding was carried out using and digestive gland tissues was also confirmed by RT-qPCR as the Native Barcoding Expansion kit (EXP-NBD104), according described above. to manufacturer’s instruction. Here, the ‘HAV + oyster’ sample was barcoded with the NB01 barcode and the ‘oyster’ sample with NB02. Barcoding was carried out in a 50 µl reaction with Direct cDNA Library Preparation and 22.5 µl of end-prep cDNA, 2.5 µl of native barcode (NB01 or Nanopore Sequencing NB02) and 25 µl of Blunt/TA ligase Master Mix (New England Library preparation was carried out using a direct Biolabs). After mixing by flicking the tube, the reaction mix was cDNA native barcoding protocol (SQK-DCS 109 with incubated for 10 min at room temperature. The barcoded cDNA EXP-NBD104) as recommended by ONT (Version: was then purified using AMPure XP beads as described in the DCB_9091_v109_revH_14Aug2019). The amount of RNA ONT protocol and resuspended in 26 µl of nuclease-free water. recommended by ONT for this protocol is 100 ng of The barcoded samples were pooled, and the volume adjusted to polyA+RNA. Two samples were prepared using: 1) a mixture 65 µl to which 5 µl of Adapter Mix II (AMII, ONT), 20 µl of 5X of RNA extracted from FRhK cells infected with HAV and NEBNext Quick Ligation Reaction Buffer (New England Biolabs) RNA extracted from oyster tissue (HAV + oyster); and 2) RNA and 10 µl of Quick T4 DNA ligase (New England Biolabs) were extracted from oyster tissues (oyster). For the preparation of the added. The reaction mix was incubated at room temperature first library, 49 ng of RNA extracted from FRhK cells infected for 10 min and purified using AMPure XP beads as described with HAV was mixed with 51 ng of oyster tissue RNA. For the in the ONT protocol and resuspended in 13 µl of nuclease-free Frontiers in Microbiology | www.frontiersin.org 3 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing water. 12 µl of the library solution were mixed with 37.5 µl of and end-prepped by mixing 0.5 µl DNA CS (ONT), 23.5 µl of sequencing buffer (SQB, ONT) and 25.5 µl of Loading Beads amplicons (ca. 140 fmol), 1.75 µl NEBNext FFPE DNA Repair (LB, ONT) and loaded into a flow cell (FLO-MIN106, ONT) Buffer (New England Biolabs), 1 µl NEBNext FFPE DNA Repair equipped with R9.4.1 chemistry on a MinION (Mk1B, ONT) Mix (New Englands Biolabs), 1.75 µl Ultra II End-prep reaction device. The MinION was operated using MinKNOW and the buffer (New England Biolabs) and 1.5 µl Ultra II End-prep flow cell was primed following manufacturer’s instructions, with enzyme mix (New England Biolabs). The mixture was mixed by a total run time of 26 h. gently flicking the tube, spin down and incubated at 20 ◦ C for 5 min and 65◦ C for 5 min. The amplicons were cleaned-up using AMPure XP beads as described in the ONT protocol and eluted Amplicon Library Preparation and in 31 µl of nuclease-free water. Adapters were then ligated by Nanopore Sequencing mixing 30 µl of amplicons prepared in the previous step, 12.5 µl For the amplicon library preparation, cDNA was synthetized of ligation buffer (ONT), 5 µl of NEBNext Quick T4 DNA ligase using RNA extracted from HAV FRhK infected cells as described (New England Biolabs) and 2.5 µl of adapter mix (ONT). The above. For the reverse transcription, 10 µl of RNA, 1 µl of mixture was mixed by gently flicking the tube, spin down and random hexamers (50 µM, Sigma), 1 µl of dNTPs (10 mM, New incubated at room temperature for 10 min. The amplicons were England Biolabs) and 1 µl of nuclease-free water were mixed and cleaned-up using AMPure XP beads and using 125 µl of fragment incubated at 65 ◦ C for 5 min and immediately after placed on ice buffer (ONT) as described in the ONT protocol and the pellet for 4 min. In a separate tube, a mixture was prepared containing was resuspended in 8 µl of elution buffer (ONT). Three µl of this 4 µl of SuperScript IV RT buffer (Invitrogen, Life Technologies), suspension were used for quantification using the Qubit DNA HS 1 µl of DTT (1 mM, Invitrogen, Life Technologies), 1 µl of Assay (Invitrogen, Life Technologies). The sequencing mix was RNaseOUT (40 U/µl, Invitrogen, Life Technologies) and 1 µl of prepared using 10 fmol of each prepared sample (VP1 or VP3) SuperScript IV (200 U/µl, Invitrogen, Life Technologies). This and eluted in elution buffer (ONT) in a total volume of 5 µl which mixture was added to the RNA-primer mix and incubated at was then mixed with 10 µl of loading beads (ONT) and 15 µl of 23◦ C for 10 min followed by 30 min at 50◦ C and inactivation sequencing buffer (ONT). The Flongle flow cell (which had 71 of the reaction at 80◦ C for 10 min. A 455-bp fragment of active nanopores) was then primed with a mix of 117 µl of flush VP3 region (corresponding to nucleotides 1470 to 1839 of HAV buffer with 3 µl of flush tether, followed by the addition of 30 µl of M59808) was amplified by PCR using the forward primer VP3- the sequencing mix. The Flongle was operated using MinKNOW 1431B (50 -CACT CAA TGT TTT AGC TAG A-30 ) and reverse and the total run time was 16 h. primer NH2 -VP3 (50 - TCT ACC TGA ATG ATA TTT GG - 30 ) described by Sanchez et al. (2003). Additionally, a 498 bp fragment of the VP1 region (corresponding to nucleotides 2394 to 2891 of HAV M59808) was amplified by PCR using the forward VP1F primer (50 - TCC TGA ATT GAA ACC TGG AGA A - 30 ) and the reverse VP1R_primer (50 - GCA ATC TGA ATG AAA CCA ATC CA -30 ), designed using primer3 software with default parameters and the consensus obtained in the present study (by directed cDNA sequencing) as the source sequence. For both VP1 and VP3 fragments amplification, PCR reactions were performed in a final reaction mix volume of 25 µl and contained 5 µl of 5X Colorless GoTaq Flexi Buffer (Promega), 12.3 µl of nuclease free water, 1 µl of dNTP (10 mM, New England Biolabs), 2 µl of MgCl2 (25 mM, Promega), 1.25 µl of each primer (10 µM), 0.2 µl of GoTaq G2 Hot Start (Promega) and 2.0 µl of template cDNA. The amplification cycle consisted of an initial denaturation step at 95◦ C for 2 min, followed by 38 cycles of 30 s at 95◦ C, 30 s at 50◦ C, 1 min at 72◦ C and with a final extension of 10 min at 72◦ C. All PCR reactions were performed in an Eppendorf Master Cycler Nexus X2, (Eppendorf). Amplicons were separated by electrophoresis on 1.5% agarose gels, stained with Ethidium bromide and visualized under UV light to confirm the expected size and absence of non-specific amplification (Figure 1). The quantity of the amplicons obtained was determined using the Qubit DNA HS Assay (Invitrogen, Life FIGURE 1 | Amplicons obtained with the primers VP1F and VP1R that amplify Technologies). Library preparation was carried out using an a 498-bp fragment of the VP1 region of HAV (Lanes 1 and 2), and amplicons obtained with the primers NH2 and VP3-1431B that amplify a 455-bp Amplicon by ligation protocol (SQK-LSK109) recommended by fragment of the VP3 region of HAV (Lanes 3 and 4). Lanes N, negative ONT (Version: DGE_9063_v109_revS_14Aug2019) but using a controls. Lanes L, 100 bp DNA ladder (Promega). Flongle device (flow cell and adapter). Amplicons were repaired Frontiers in Microbiology | www.frontiersin.org 4 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing Bioinformatic Analysis to the low concentration obtained in our samples (7 ng/µL, Base calling of Fast5 files was performed in a MinIT (ONT- Qubit quantification). For oyster RNA (1562 ng/µL, Qubit minit-release 19.05.2) with Flip Flop algorithm. Demultiplexing quantification), the nanodrop OD 260/280 and 260/230 ratios and removal of barcode sequences was performed using were 2.11 and 2.19, respectively. These values are indicative Porechop (v0.2.3, https://github.com/rrwick/Porechop). Since a of pure RNA with no indication of contaminants. We have large proportion of the reads were not successfully barcoded, tested other RNA extraction and purification methods (data not all reads were pooled for further analysis. Run metrics were shown) but these have often resulted in lower 260/230 ratios calculated using the Nanoplot (v1.20.0) and reads were filtered comparatively to the TRIzol/Zymo kit used. Other methods were (sequence quality > 7) using Nanofilt (v2.3.0) (De Coster et al., therefore not used for RNA preparation in this work. The HAV 2018). Filtered fastq files were converted to fasta files, which level quantified by RT-qPCR was 1.7 × 108 copies/mL of cell- were then used to create a cDNA nucleotide database using culture supernatant fluid. The recommended quantity of RNA Blast v2.9.0+ (Camacho et al., 2009). We used as query a for the direct cDNA ONT protocol (SQK-DCS109) is 100 ng hepatitis A virus genome (strain HM-175, wild type, Genbank of RNA with a polyA tail. As the total RNA extracted from accession number M14707) which was blasted against our FRhK cells infected with HAV for library preparation was below cDNA nucleotide database to detect HAV reads. Reads showing this value (i.e., a total of 49 ng obtained from 2.5 ml of viral a significant e-value were extracted from the original fasta suspension) an additional 51 ng of oyster RNA was added. files and screened for chimeras using YACRD, which resulted Based on the RT-qPCR results, we estimated that ca. 4 × 108 in 302 HAV chimera-free reads (Marijon et al., 2019). The viral genome copies were used for the library preparation which extracted HAV reads were mapped against the HAV genome should equate to approximately 1.81 ng of viral RNA (0.71 fmol) (Genbank accession number M14707) using Minimap2 (Li, 2018) assuming a genome size of 7.5 kilobases. Therefore, the estimated and the alignment was visualized using Tablet (Milne et al., amount of eukaryotic RNA extracted from the FRhK cells and 2013). The HAV chimera-free 302 reads were assembled using oyster gills was 288.4 fmol and hence 406 times more than Canu (v1.8) using default parameters for nanopore sequence viral RNA copies. Most importantly, such limited amount of data (Koren et al., 2017) to produce a consensus nucleotide viral RNA often represents a limiting factor for RNA or cDNA sequence. Moreover, a consensus sequence was obtained by direct sequencing approaches (Posada-Cespedes et al., 2017). mapping the HAV reads to HAV cytopathogenic pHM175/18f In clinical and food samples PCR-based methods are currently strain genome (accession number M59808) using LAST sequence the only possible approach to sequence viral RNA (Acevedo aligner (version 923, with default parameters) as implemented and Andino, 2014). Regardless, whole-genome sequencing using in Nanopipe (Shabardina et al., 2019). Multiple alignments were PCR-based approaches can also be quite laborious and are carried out using MAFFT v7 software (Katoh and Standley, known to induce several PCR amplification artifacts and bias 2013) using default parameters, with the nucleotide sequences (Ardui et al., 2018). of the consensus obtained in the present study using Canu and the genomes of HAV genotypes IA (KC182590, EF406357, Direct cDNA Nanopore Sequencing X75215, X83303), IB (AF314208, M20273, M14707, AF268396, For the direct cDNA approach, a total of 2,792,673 nanopore raw M59808), IIA (AY644676), IIB(AY644670), IIIA (AJ299464), IIIB reads were obtained with a N50 (i.e., defined as the minimum (AB258387) as well as the simian HAV genome (EU140838). The contig length needed to cover 50% of the genome) of 1,044 bp sequence alignment region used for downstream analysis, was and mean read length of 946 bp. The mean read quality of retained from position 84 up to 7422 regions of the HAV genome unfiltered reads was 7.5 and approximately 67% of the reads had M14707. Maximum Likelihood (ML) analysis was conducted on the sequence alignment region shared by all genome sequences using MEGA7 (Kumar et al., 2016) using the General Time Reversible model (Nei and Kumar, 2000), gamma distribution with invariant sites (G+I), four discrete rate categories and 1,000 bootstrap replicates. Single nucleotide variants (SNVs) were called using Longshot software with a strand bias p-value of 0.01 (-c 50, -C 100,000) (Edge and Bansal, 2019). RESULTS AND DISCUSSION RNA Extraction and Purification An important factor that impacts cDNA synthesis from RNA viruses and ONT library preparation is the quality of the RNA and the absence of contaminants that can be carried FIGURE 2 | Alignment length and nucleotide sequence identity of raw over during the RNA extraction. In our experiments, these nanopore hepatitis A virus reads identified by BLASTn analysis using the parameters could not be fully explored using the nanodrop OD genome of HAV wild type (M14707) as a reference. 260/280 and 260/230 ratios for the FRhK cells/viral RNA due Frontiers in Microbiology | www.frontiersin.org 5 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing a quality above 7 (Q7). After eliminating reads with a quality of reads was only 0.013% which is considerably less than the below 7 a total of 1,860,348 reads were retained. A BLASTn expected 0.245% based on the estimated quantification by RT- search using the HAV genome (M14707) as the query and the qPCR described above. The current protocol also required the Q7 reads as the database allowed the identification of 243 HAV presence of a poly A-tailed RNA and thus cDNA will only be reads (before chimera analysis) with a significant e-value (<1E- synthesized for RNA containing a poly A tail, which in HAV it 23). The percentage of HAV reads relative to the total number stretches at the 30 end (Daijogo and Semler, 2011). The relatively FIGURE 3 | Alignment of HAV nanopore reads obtained by direct cDNA sequencing using the genome of HAV wild type (M14707) as a reference. Detail of part of the alignment of 16 HAV nanopore reads highlighting the deletions (marked with a red asterisk) and substitutions. FIGURE 4 | Phylogenomic analysis using HAV genomes (several Human genotypes IA, IB, IIA, IIB, IIIA, and IIIB as well as a Simian HAV) and the consensus sequence obtained in the present study (marked with a red arrow). Bootstrap values obtained by Maximum Likelihood analysis are indicated at the main nodes. Frontiers in Microbiology | www.frontiersin.org 6 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing FIGURE 5 | (A) Position and frequency of nucleotides of candidate sequence nucleotide variants (SNVs) detected using the HAV reads obtained by ONT direct cDNA sequencing. (B) Frequency of nucleotides observed in the VP1 amplicon dataset. The position of the nucleotides is shown using the HAV cytopathogenic pHM175/18f strain genome as reference (accession number M59808). low number of HAV reads detected can be partially explained the assembly was performed using the HAV strain pHM175/18f by RNA shearing during extraction and library preparation. (M59808) as reference a consensus with a sequence similarity of Notwithstanding, several long sequence reads (with more than 99.8 % was obtained (7423 bases, Genbank accession number 7Kb) representing HAV sequences were identified, suggesting MT181522). The alignment of this reference-based consensus that some viral particles were still intact. The sequence identity with HAV M59808 strain revealed 12 substitutions and 1 of the HAV reads ranged between 81.1 and 96.6% with a ambiguous position. Since no other sequencing methods (e.g., mean value of 88.9% and a standard deviation of 2.8%, when Illumina sequencing, which has a much lower error rate) were using the HAV wild type genome (M14707) as a reference used in this study, it was not possible to determine the sequencing (Figure 2). The length of these reads ranged between 129 and errors accurately. Nevertheless, the 0.2% nucleotide sequence 7667 bases with a median value of 1173 bases. When analyzing dissimilarity observed between the consensus and the HAV only the long reads that had more than 7.3 kilobases, these cytopathogenic pHM175/18f strain suggests that it is possible contained almost the whole genome of the HAV. In fact, when to circumvent ONT sequencing errors by using fine-tuned data comparing the 7.3 Kb reads with the HAV reference wild type analysis pipelines. Notwithstanding, RNA viruses exhibit high (M14707) the nucleotide sequence identity was between 84.7 levels of mutation rates (Duffy et al., 2008) and this can also and 90.6% with indels and substitutions accounting for a mean contribute to the presence of new viral strains during replication value of 5.6 and 5.9% of the dissimilarity, respectively. The via point mutations, insertions and deletions (Posada-Cespedes resulting reference-based mapping covered the whole genome et al., 2017), especially in strains used in laboratories after many of the HAV (Figure 3). One of the major limitations of passages and cultivation. nanopore sequencing is still its relatively high error rates, with the proportion of read errors often ranging from 10-30% (Goodwin et al., 2015). Nonetheless several algorithms have HAV Genetic Diversity been developed to reduce these sequencing errors (Goodwin The high error rate of nanopore sequencing makes the detection et al., 2015; Salmela et al., 2017). In order to circumvent of single nucleotide variants quite challenging (Vasudevan et al., reference-based read alignment biases, we then utilized the 2019). These technical limitations present in Oxford nanopore Canu software that operates by improving the accuracy of sequence data can be bypassed with bioinformatic tools that bases in reads, trim reads to the portion that have high- enable to distinguish artifacts from true polymorphisms, such quality sequence followed by de novo assembly. From this is the case of the newly developed Longshot software (Edge de novo assembly a consensus nucleotide sequence with 7362 and Bansal, 2019). In the present study, 7 possible SNVs bases was obtained (after trimming the poly A tail and the were detected using Longshot software being 6 transitions and 30 -end using the M14707 HAV genome as reference). This 1 transversion (Figure 5). Previous studies have also reported consensus showed a nucleotide sequence identity of 99.0% (7349 a considerably higher number of transitions over transversions out of 7423 identical bases) against the genome of a HAV in HAV (Sanchez et al., 2003). These results suggest that the cytopathogenic pHM175/18f strain (accession number M59808). HAV population characterized might have been composed by We identified 61 deletions and 13 substitutions along the 7423 two or more haplotypes. This highlights the potential of direct bases alignment of the consensus with the HAV cytopathogenic cDNA nanopore sequencing to characterize the genetic diversity pHM175/18f strain (M59808). Phylogenetic analysis of the HAV of RNA viruses, which can be an alternative to PCR based consensus obtained using Canu and reference sequences of HTS methods and its inherent amplification biases (Posada- known genotypes showed that the consensus clustered robustly Cespedes et al., 2017). To confirm the presence of SNVs detected with the sub-genotype IB group (Figure 4). However, when using direct cDNA sequencing, we used an ONT amplicon-based Frontiers in Microbiology | www.frontiersin.org 7 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing approach targeting two regions of the HAV genome (VP1 and sequencing approach and obtained long-reads some of which VP3). For the VP3 amplicon dataset, a total of 192, 738 reads were over 7.3 kilobases in length, containing nearly the whole were mapped against HAV M59808 strain. The 452-bp consensus genome of the virus. sequence obtained (without the annealing regions) showed 100% nucleotide sequence similarity with the consensus obtained using the cDNA approach. Moreover, no polymorphisms were DATA AVAILABILITY STATEMENT detected in the VP3 amplicon dataset which is in accordance with the results obtained using the direct cDNA dataset. For The datasets generated for this study can be found in Genbank, the VP1 amplicon dataset, a total of 124,986 reads were mapped under accession number MT181522. against the HAV M59808 strain. One hundred percent nucleotide sequence similarity was also observed between the 415-bp consensus sequence obtained using the amplicon data and the AUTHOR CONTRIBUTIONS consensus obtained using the direct cDNA sequencing approach. However, no polymorphisms were detected in the VP1 amplicon FB and JM-U conceived and designed the study. FB, TS, JL, dataset which contrasts with the 3 SNVs detected in the same RS, and CP did the laboratory work. FB, VF, and RvA did the region using the direct cDNA sequencing approach (positions bioinformatic analyses. FB, TS, JL, VF, CP, DW, RvA, and JM-U 2596, 2685, and 2751, Figure 5). This could be due to the wrote the manuscript. fact that amplification steps can exponentially amplify the most abundant viral clones, resulting in scarce or absent representation of low abundance viral haplotypes (Posada-Cespedes et al., 2017). FUNDING In this regard, our PCR-free approach could circumvent these amplification bias. This study was funded by Cefas Seedcorn programme (references 40200FT and DP422). CONCLUSION ACKNOWLEDGMENTS In the present study we demonstrated the potential of nanopore single-molecule sequencing to sequence the whole genome We would like to thank the staff of the Food Safety Group for of HAV from infected cell lines. We were able to sequence providing laboratory support. We would also like to thank two HAV by direct cDNA sequencing using a PCR-free nanopore reviewers for their suggestions and comments. REFERENCES Edge, P., and Bansal, V. (2019). Longshot enables accurate variant calling in diploid genomes from single-molecule long read sequencing. Nat Commun 10, 4660. Acevedo, A., and Andino, R. (2014). Library preparation for highly accurate doi: 10.1038/s41467-019-12493-y population sequencing of RNA viruses. Nat Protoc 9, 1760–1769. doi: 10.1038/ Goodwin, S., Gurtowski, J., Ethe-Sayers, S., Deshpande, P., Schatz, M. C., nprot.2014.118 and McCombie, W. R. (2015). Oxford Nanopore sequencing, hybrid error Andino, R., and Domingo, E. (2015). Viral quasispecies. Virology 47, 46–51. doi: correction, and de novo assembly of a eukaryotic genome. Genome Res 25, 10.1016/j.virol.2015.03.022 1750–1756. doi: 10.1101/gr.191395.115 Ardui, S., Ameur, A., Vermeesch, J. R., and Hestand, M. S. (2018). Single molecule Hansen, S., Dill, V., Shalaby, M. A., Eschbaumer, M., Bohlken-Fascher, S., real-time (SMRT) sequencing comes of age: applications and utilities for Hoffmann, B., et al. (2019). Serotyping of foot-and-mouth disease virus using medical diagnostics. Nucleic Acids Research 46, 2159–2168. doi: 10.1093/nar/ oxford nanopore sequencing. J Virol Methods 263, 50–53. doi: 10.1016/j. gky066 jviromet.2018.10.020 Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. Katoh, K., and Standley, D. M. (2013). MAFFT multiple sequence alignment (2009). BLAST+: architecture and applications. BMC Bioinformatics 10:421. software version 7: improvements in performance and usability. Mol Biol Evol doi: 10.1186/1471-2105-10-421 30, 772–780. doi: 10.1093/molbev/mst010 Costafreda, M. I., Bosch, A., and Pinto, R. M. (2006). Development, evaluation, Keller, M. W., Rambo-Martin, B. L., Wilson, M. M., Ridenour, C. A., Shepard, and standardization of a real-time TaqMan reverse transcription-PCR S. S., Stark, T. J., et al. (2018). Direct RNA Sequencing of the Coding Complete assay for quantification of hepatitis A virus in clinical and shellfish Influenza A Virus Genome. Scientific Reports 8, 14408. doi: 10.1038/s41598- samples. Appl Environ Microbiol 72, 3846–3855. doi: 10.1128/AEM.02 018-32615-8 660-05 Koren, S., Walenz, B. P., Berlin, K., Miller, J. R., Bergman, N. H., and Phillippy, Daijogo, S., and Semler, B. L. (2011). “Chapter 1 - Mechanistic Intersections A. M. (2017). Canu: scalable and accurate long-read assembly via adaptive Between Picornavirus Translation and RNA Replication,” in Advances in k-mer weighting and repeat separation. Genome Res 27, 722–736. doi: 10.1101/ Virus Research, eds K. Maramorosch, A. J. Shatkin, and F. A. Murphy gr.215087.116 (Cambridge, MA: Academic Press), 1–24. doi: 10.1016/B978-0-12-385987-7. Kumar, S., Stecher, G., and Tamura, K. (2016). MEGA7: Molecular Evolutionary 00001-4 Genetics Analysis Version 7.0 for Bigger Datasets. Mol Biol Evol 33, 1870–1874. De Coster, W., D’Hert, S., Schultz, D. T., Cruts, M., and Van Broeckhoven, doi: 10.1093/molbev/msw054 C. (2018). NanoPack: visualizing and processing long-read sequencing data. Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 2666–2669. doi: 10.1093/bioinformatics/bty149 Bioinformatics 34, 3094–3100. doi: 10.1093/bioinformatics/bty191 Duffy, S., Shackelton, L. A., and Holmes, E. C. (2008). Rates of evolutionary change Marijon, P., Chikhi, R., and Varré, J.-S. (2019). Yacrd and fpa: upstream in viruses: patterns and determinants. Nat Rev Genet 9, 267–276. doi: 10.1038/ tools for long-read genome assembly. bioRxiv[Preprint] doi: 10.1101/67 nrg2323 4036 Frontiers in Microbiology | www.frontiersin.org 8 May 2020 | Volume 11 | Article 874 Batista et al. HAV Nanopore Sequencing McNaughton, A. L., Roberts, H. E., Bonsall, D., de Cesare, M., Mokaya, J., Somerville, V., Lutz, S., Schmid, M., Frei, D., Moser, A., Irmler, S., et al. (2019). Lumley, S. F., et al. (2019). Illumina and Nanopore methods for whole genome Long-read based de novo assembly of low-complexity metagenome samples sequencing of hepatitis B virus (HBV). Sci Rep 9, 7081. doi: 10.1038/s41598- results in finished genomes and reveals insights into strain diversity and 019-43524-9 an active phage system. BMC Microbiology 19:143. doi: 10.1186/s12866-019- Milne, I., Stephen, G., Bayer, M., Cock, P. J., Pritchard, L., Cardle, L., et al. (2013). 1500-0 Using Tablet for visual exploration of second-generation sequencing data. Brief Takeda, H., Yamashita, T., Ueda, Y., and Sekine, A. (2019). Exploring the hepatitis C Bioinform 14, 193–202. doi: 10.1093/bib/bbs012 virus genome using single molecule real-time sequencing. World J Gastroenterol Nei, M., and Kumar, S. (2000). Molecular Evolution and Phylogenetics. New York, 25, 4661–4672. doi: 10.3748/wjg.v25.i32.4661 NY: Oxford University Press. Tan, C. C. S., Maurer-Stroh, S., Wan, Y., Sessions, O. M., and de Sessions, P. F. Ng, S. H., Braxton, C., Eloit, M., Feng, S. F., Fragnoud, R., Mallet, L., et al. (2018). (2019). A novel method for the capture-based purification of whole viral native Current Perspectives on High-Throughput Sequencing (HTS) for Adventitious RNA genomes. AMB Express 9, 45. doi: 10.1186/s13568-019-0772-y Virus Detection: Upstream Sample Processing and Library Preparation. Viruses Treangen, T. J., and Salzberg, S. L. (2011). Repetitive DNA and next-generation 10, 566. doi: 10.3390/v10100566 sequencing: computational challenges and solutions. Nature reviews. Genetics Posada-Cespedes, S., Seifert, D., and Beerenwinkel, N. (2017). Recent advances in 13, 36–46. doi: 10.1038/nrg3117 inferring viral diversity from high-throughput sequencing data. Virus Res 239, Ueda, Y. (2019). New sequencing technology, third-generation Nanopore 17–32. doi: 10.1016/j.virusres.2016.09.016 sequencing, for hepatitis C virus research. Hepatol Res. 49, 961–963. doi: 10. Quinones-Mateu, M. E., Avila, S., Reyes-Teran, G., and Martinez, M. A. (2014). 1111/hepr.13383 Deep sequencing: becoming a critical tool in clinical virology. J Clin Virol 61, Vasudevan, K., Devanga Ragupathi, N. K., Jacob, J. J., and Veeraraghavan, B. (2019). 9–19. doi: 10.1016/j.jcv.2014.06.013 Highly accurate-single chromosomal complete genomes using IonTorrent and Salmela, L., Walve, R., Rivals, E., and Ukkonen, E. (2017). Accurate self-correction MinION sequencing of clinical pathogens. Genomics 112, 545–551. doi: 10. of errors in long reads using de Bruijn graphs. Bioinformatics 33, 799–806. 1016/j.ygeno.2019.04.006 doi: 10.1093/bioinformatics/btw321 WHO (2017). Global Hepatitis Report. Geneva: World Health Organization. Sanchez, G., Bosch, A., Gomez-Mariano, G., Domingo, E., and Pinto, R. M. (2003). Evidence for quasispecies distributions in the human hepatitis A virus genome. Conflict of Interest: The authors declare that the research was conducted in the Virology 315, 34–42. doi: 10.1016/s0042-6822(03)00483-5 absence of any commercial or financial relationships that could be construed as a Sattar, S. A., Jason, T., Bidawid, S., and Farber, J. (2000). Foodborne spread of potential conflict of interest. hepatitis A: Recent studies on virus survival, transfer and inactivation. The Canadian journal of infectious diseases 11, 159–163. doi: 10.1155/2000/805156 Copyright © 2020 Batista, Stapleton, Lowther, Fonseca, Shaw, Pond, Walker, van Schatz, M. C., Delcher, A. L., and Salzberg, S. L. (2010). Assembly of large genomes Aerle and Martinez-Urtaza. This is an open-access article distributed under the using second-generation sequencing. Genome Research 20, 1165–1173. doi: 10. terms of the Creative Commons Attribution License (CC BY). The use, distribution 1101/gr.101360.109 or reproduction in other forums is permitted, provided the original author(s) and Shabardina, V., Kischka, T., Manske, F., Grundmann, N., Frith, M. C., Suzuki, Y., the copyright owner(s) are credited and that the original publication in this journal et al. (2019). NanoPipe-a web server for nanopore MinION sequencing data is cited, in accordance with accepted academic practice. No use, distribution or analysis. Gigascience 8, giy169. doi: 10.1093/gigascience/giy169 reproduction is permitted which does not comply with these terms. Frontiers in Microbiology | www.frontiersin.org 9 May 2020 | Volume 11 | Article 874