Molecular characterization of varicella zoster virus isolated from clinical samples in India

Varicella zoster virus (VZV) is an alpha-herpesvirus and a causative agent of chickenpox (varicella) and shingles (zoster). Chickenpox is the primary infection characterized by fever and a generalized rash, which is prevalent in children. Although a live-attenuated vaccine has been available since 19741, chickenpox remains a primary childhood infection. Propagation of VZV in cell culture was first reported in 1953, and a complete nucleotide sequence for VZV was first determined using the M13-dideoxynucleotide technology from the Dumas strain2.

The genome of VZV is approximately 125 kb in length and consists of a double-stranded DNA and has 72 known open reading frames (ORFs)2. Phylogenetic analysis of VZV sequences shows that this genetically stable virus shows sufficient variation to classify it into six clades. Restriction fragment length polymorphism (RFLP) in ORF38, ORF54 and ORF62345 has been the most common method employed for virus classification. Another method classifies the virus based on scattered single-nucleotide polymorphisms (SNPs) in ORF1, ORF21, ORF50 and ORF54, thereby classifying the VZV in four main viral clades, namely A, B, C and J6. Classification of VZV by glycoprotein sequences has led to another system of classification into four significant clades designated as A, B, C and D4. The fourth system of VZV classification refers to ORF22-, ORF21- and ORF50-based genotyping and classifies the virus into five confirmed clades E1, E2, J, M1 and M2, and two provisional clades M3 and M47. Classification of VZV based on complete genome leads to six significant clades that have a distinct geographic distribution pattern89. Based on the full genome sequence-based classification, clade 2 is found to be dominantly present in Asia4. However, in countries belonging to the South Asian subcontinent, namely India, Pakistan, Bangladesh and Sri Lanka, clades 4 and 5 are found to be more prevalent56.

The genetic diversity among VZV sequences is 0.00063, which is almost 40 times less than that of other large double-stranded DNA viruses such as herpes simplex virus and cytomegalovirus8. However, with such a small genetic diversity, the clinical condition manifested by VZV is varied1011. A study from India based on clinical-epidemiological details of VZV-infected cases reported symptoms varying from mild fever with a rash to very severe (bleeding manifestation) that led to mortality10. VZV infection causes a communicable disease that mostly spreads through the vesicular skin cell, a place where it is highly concentrated12. Although VZV infections are prevalent across India1314, complete genome sequences from India are not available. The present study was thus aimed to identify the circulating clade of VZV in India using its genomic sequence from clinical samples.

Material & Methods

This study was conducted at Maximum Containment Laboratory of the Indian Council of Medical Research - National Institute of Virology (ICMR-NIV), Pune, India during April 2017. The samples used in this study were received as part of the focal VZV outbreaks that had occurred in Jaipur, Rajasthan, and Mumbai, Maharashtra, India. Written informed consent was taken before the collection of the clinical specimen from VZV-infected patients. The study was approved by the Institute’s Human Ethics Committee.

In the earlier published work, 79 samples from Maharashtra (n=50), Gujarat (n=12) and Rajasthan (n=17) States, India, were tested by real-time PCR to test the VZV positivity15. Sixty one positive clinical samples were observed from these three locations15. However, based on the availability of the sample and volume required to generate complete genomic sequences, using next-generation sequencing (NGS), two samples each from Rajasthan and Maharashtra were selected for this study.

Next-generation sequencing (NGS):

DNA library preparation: The whole-genome sequencing of VZV was performed using the NGS platform. DNA was extracted from 1 ml of the VZV-positive samples using a Qiagen DNA extraction kit (Qiagen, Hilden, Germany). The extracted DNA was quantified using Qubit 2.0 (Invitrogen, Carlsbad, USA) and diluted to 0.2 ng/μl. The diluted DNA was used to prepare the DNA library by Nextera DNA library preparation (Illumina, California, USA) using the previously described protocol16. The DNA library was quantified using the KAPA library quantification kit (Kapa Biosystems, Roche, Basel, Switzerland). The final DNA libraries to be sequenced were normalized and loaded onto the MiniSeq mid-output cartridges (300 cycles) for paired-end sequencing using the Illumina MiniSeq platform, (Illumina, USA).

The library for DNA sample with accession number MH499466 was made using NEBNext^® Ultra™ DNA Library Prep Kit for Illumina^® (New England Biolabs^® Inc. Massachusetts, USA) to retrieve more sequences. The library was made as per the manual of the kit with slight modifications. The final input concentration of the DNA sample for NGS library preparation was 100 ng. The concentration was decided based on the previous optimization performed using different DNA concentrations. Six PCR cycles were used for the amplification of adaptor-ligated DNA.

Analysis of the reads: The FastQ files generated from the Illumina machine were imported to CLC Genomics Workbench v10.1.1 (Qiagen, Aarhus, Denmark) for analysis. Initially, contiguous sequences (contigs) were generated using the de novo contig assembly method, as implemented in the CLC Genomics Workbench. These contigs were identified using the BLAST programme to be similar to the clade 5 sequences. VZV sequence belonging to clade 5 (accession number DQ457052.1) was downloaded as a reference sequence from GenBank. Further, reference-based mapping was performed to retrieve the complete genome of VZV from the FastQ files. These sequences were submitted to the public repository.

Phylogenetic analysis: Complete reference sequences for VZV genomes were downloaded from the GenBank database to generate the phylogenetic tree (Table I). These sequences were aligned and manually checked for its correctness. The difference in the nucleotide and the per cent identity was retrieved using CLC Genomics Workbench. Maximum likelihood tree was generated using GTR+ G+I as the nucleotide substitution model in MEGA version 7.0 software17. A bootstrap replication of 1000 cycles was used to assess the significance of the phylogenetic tree generated.

Results

Next-generation sequencing: Reference-based mapping using the reference genome [DQ457052.1 (CA123)] led to the retrieval of 124,700 to 121,928 bp of VZV sequences from different clinical samples (Table II). The Q30 was >95 per cent for the reads generated from the clinical samples. VZV sequences having accession numbers MH499466 and MH499467 were collected from Jaipur, Rajasthan, and MH499469 and MH499468 were collected season from Mumbai, Maharashtra. Samples with accession numbers MH499468 and MH499467 had partial genome (~97.7-98.5% genome coverage) retrieval due to the availability of lesser sample volume. The details regarding the sample type, total read generated and length of the sequences recovered along with the percentage of the reads mapped are given in Table II.

Phylogenetic analysis: Fig. 1A and B depicts the phylogenetic tree for the VZV sequences that have genomic coverage of more than 99 and 97 per cent, respectively, to reference sequence. The analysis revealed that the Indian sequences belonged to clade 5. However, an interestingly separate branch/cluster was observed with bootstrap support of above 90 and 80 per cent for the two different trees generated (Fig. 1A and B). We also determined the differences and the per cent identity of the VZV sequences analyzed in the study (Fig. 2). The upper right matrix and a lower left matrix of Fig. 2 denote the nucleotide differences and per cent identity, respectively, between the VZV reference and Indian VZV sequences obtained from Indian samples. It was observed that VZV sequences belonging to clade 5 had nucleotide per cent identity that distinctly formed two different groups, except MH499468, which might be due to its smaller length. Further, the nucleotide difference of strain CA123 (accession number: DQ457052.1) belonging to clade 5 and Indian sequences with other reference VZV sequences was also significant.

Close

Fig. 1

A maximum likelihood evolutionary tree for varicella zoster virus genomic coverage of more than 99 per cent (A) and more than 97 per cent (B) was drawn using the general time-reversible model, having a gamma parameter as 0.1000 with an invariable rate variation of 49.64 and 49.63, respectively. A bootstrap replication of the 1000 cycles was performed to assess the statistical robustness.

Export to PPT

Close

Fig. 2

The nucleotide differences and per cent identity between the varicella zoster virus sequences used in this study obtained from the CLC Genomics Workbench. Nucleotide difference is shown in the upper right panel and per cent identity on the left lower panel.

Export to PPT

Table III enlists the 13 SNPs observed in the Indian VZV five sub-clade compared to reference clade 5 genomes (accession no.: DQ457052.1). Three of the 13 SNPs (A72484C, G74974A and C75130T) and a single multi-nucleotide polymorphism (MNP) (CC75126-75127TT) were observed in the ORF40 region that encodes for major capsid protein. Of the remaining nine SNPs, non-synonymous amino acid changes were observed for four SNPs (Table III).

Discussion

The sequencing of large genomes like VZV is a challenging task. In this study the NGS platform was used to sequence the genome as (i) VZV being a slow-growing virus in tissue culture, the isolations took long time, and (ii) clinical samples were limited. Despite these, nearly complete genomes of VZV were obtained using NGS, albeit with some gaps in the sequences. All the four VZV sequences formed a cluster with clade 5 sequences. The Indian VZV sequences formed a distinct cluster with a separate branch having 70 per cent bootstrap support. This indicated that a different sub-clade belonging to clade 5 was circulating in the Rajasthan and Maharashtra States of India. As the sample size was restricted, the existence of a clade other than the observed sub-clade could not be denied.

Depledge et al18 in their study tried to understand if there was a correlation between viral genetic diversity and disease progression in terms of central nervous system (CNS) and non-CNS cases. They observed that cases with encephalitis showed more considerable genetic variability as compared to VZV-infected cases with meningitis, usually a self-limiting infection. In another study by Peters et al19, it was hypothesized that the virus must have co-evolved along with its host, and the recombination among various strains with the host genome might be responsible for variations specific to a given geographic region. All four genome sequences obtained from India belonged to clade 5, and formed a separate branch, thereby suggesting that there might have been geographical region-specific host interactions leading to various nucleotide changes.

Sequence data from India are fragmented. The present analysis helped to understand the phylogenetic analysis of nearly 98 per cent complete VZV sequences from India and its comparison with other global VZV sequences. The changes in the structure or function of the major capsid protein caused from the change in the SNPs needs to be further experimentally studied. Attempt to correlate disease symptoms and specific sequence variations in the VZV genome was none, apart the earlier published work18. The SNPs observed in this study were from the limited dataset, which was also the limitation of this study, and indicated a need for VZV sequencing from a substantial number of confirmed clinical cases. This will help in understanding the viral diversity associated with the clinical manifestation in the disease condition.