Lamivudine/Adefovir Treatment Increases the Rate of Spontaneous Mutation of Hepatitis B Virus in Patients

The high levels of genetic diversity shown by hepatitis B virus (HBV) are commonly attributed to the low fidelity of its polymerase. However, the rate of spontaneous mutation of human HBV in vivo is currently unknown. Here, based on the evolutionary principle that the population frequency of lethal mutations equals the rate at which they are produced, we have estimated the mutation rate of HBV in vivo by scoring premature stop codons in 621 publicly available, full-length, molecular clone sequences derived from patients. This yielded an estimate of 8.7 × 10−5 spontaneous mutations per nucleotide per cell infection in untreated patients, which should be taken as an upper limit estimate because PCR errors and/or lack of effective lethality may inflate observed mutation frequencies. We found that, in patients undergoing lamivudine/adefovir treatment, the HBV mutation rate was elevated by more than sixfold, revealing a mutagenic effect of this treatment. Genome-wide analysis of single-nucleotide polymorphisms indicated that lamivudine/adefovir treatment increases the fraction of A/T-to-G/C base substitutions, consistent with recent work showing similar effects of lamivudine in cellular DNA. Based on these data, the rate at which HBV produces new genetic variants in treated patients is similar to or even higher than in RNA viruses.


Introduction
With around 250 million people infected worldwide, hepatitis B virus (HBV) constitutes a major cause of cirrhosis, liver failure and hepatocellular carcinoma [1]. HBV has a small (3.2 kb), partially double-stranded, circular DNA virus that is primarily encapsidated as a pre-genomic RNA that undergoes reverse transcription after encapsidation [2,3]. The HBV genome is extremely compact, with four, partially or fully overlapped open reading frames (C, S, P, and X). The preC and C reading frames encode the e antigen (HBeAg) and core protein (HBcAg) respectively, the preS and S reading frames encode three forms of surface proteins (HBsAg) sharing the C-terminus (small, middle and large), the P gene encodes the viral polymerase, which acts both as reverse transcriptase (RT) and DNA-dependent DNA polymerase, and the X gene encodes a transcriptional trans-activator protein. HBV treatment consists of immune modulators (interferon) combined with nucleoside analogues whose primary effect is to inhibit the viral polymerase [4,5]. However, treatment typically fails to clear the virus, which remains active in the form covalently closed circular DNA within hepatocyte nuclei [6], and long-term treatment with nucleoside analogues often selects for drug-resistant mutations in the HBV polymerase [7].
HBV shows a molecular evolutionary rate on the order of 10 −4 nucleotide substitutions per site per year (s/s/y) at the epidemiological scale [8,9], and of 10 −5 -10 −4 s/s/y at the intra-patient level [10,11]. This rate is notably higher than those of most non-reverse transcribing DNA viruses, and similar to those of RNA viruses and retroviruses [12]. It is commonly accepted that such fast evolution is determined ultimately by the low replication fidelity the HBV polymerase which, similar to other RTs, lacks 3´exonclease proofreading activity [13] and produces frequent replication errors, resulting in highly diverse viral populations [14][15][16][17]. Despite growing interest in studying the relationship between HBV genetic variation and viral pathogenesis, as well as its implications for viral detection, prevention, treatment and prognosis, the rate of spontaneous mutation of human HBV has not been determined. Notice that, whereas the molecular evolutionary rate refers to a population genetics process, the rate of spontaneous mutation refers to a biochemical process and, therefore, these two rates should not be confused. The former describes the fixation of new alleles in a population and is determined by factors such as selection and random drift, whereas the latter is defined as the probability that new mutations appear per round of genome copying, or per cell infection cycle, and is determined by factors such as polymerase fidelity, DNA/RNA editing, and spontaneous damage. Early work with the related duck HBV estimated the reversion frequency of a single deleterious Gto-A nucleotide substitution [18] and, after correcting for several possible confounders such selection and the unknown number of rounds of replication, this yielded a mutation rate estimate of 2 × 10 −5 per nucleotide per cell infection cycle (m/n/c) [19]. However, the reliability of this value was compromised by the fact that only a single genome site was analyzed.
Another unaddressed question is the effect of treatment on the HBV mutation rate. Broadspectrum nucleoside analogues such as ribavirin are known to mutagenize RNA viruses [20][21][22][23][24][25]. Similarly, nucleoside analogue RT inhibitors such as lamivudine and AZT have a direct mutagenic effect on HIV-1 in cell culture, in addition to their inhibitory effects [26]. Moreover, treatment with nucleoside analogues can select for resistance mutations that modify the replication fidelity of the HIV-1 RT, producing an additional, indirect effect on the viral mutation rate [27,28]. In HBV, analysis of nucleotide misincorporation kinetics in vitro showed that typical lamivudine-resistant polymerase variants such as M204I and M204V display increased replication fidelity in the absence of the drug [13], suggesting a mutagenic effect for lamivudine in patients. Nucleoside analogue RT inhibitors including lamivudine have also been found to be mutagenic in different cell lines, as shown by analysis of reporter cellular genes [29,30].
Rates of spontaneous mutation can be inferred in vivo from intra-host sequence diversity data using the lethal mutation method, which is based on the principle that the population frequency of mutations that abolish viral infectivity (lethal mutations) equals the rate at which they are produced, as these cannot be inherited [31]. Using premature stop codons in essential genes as a proxy to lethal mutations, we and others have previously used this method for estimating the mutation rate of hepatitis C virus (HCV) [21,32] and HIV-1 [33] from patientderived sequences. Here, we applied the lethal mutation method to publicly available, fulllength, molecular clone HBV sequences derived from untreated patients, as well as from patients undergoing long-term lamivudine/adefovir combination therapy. We obtained an estimated spontaneous rate on the order of 10 −5 m/n/c for HBV, albeit with ample variations among HBV reading frames. Furthermore, we show that lamivudine/adefovir elevates the spontaneous HBV mutation rate by at least sixfold. Analysis of intra-host sequence diversity indicated that lamivudine/adefovir treatment produces a shift in the HBV mutation spectrum whereby A-to-G and the complementary T-to-C base transitions are specifically elevated.

Material and Methods Data
We searched molecular clone full-length sequences deposited in GenBank. Publications were first browsed in PubMed based on the following electronic search strategy: ((hepatitis b virus) AND genome Ã ) AND (((clone Ã ) OR cloning) OR full-length). This allowed us to retrieve a large number of articles, which were then inspected to collect only sequences that explicitly included the method used to obtain HBV molecular clones. Also, if fewer than three sequences per patient were analyzed, these were discarded. All sequences corresponded to plasma or serum samples, except for AF182805 and AF182804, which were obtained from a resected peritumor liver tissue and a hepatocellular carcinoma, respectively. Sequences were aligned using the MUSCLE algorithm (www.drive5.com/muscle) and refined further by visual inspection. A neighbor-joining tree was obtained with MEGA6 (www.megasoftware.net) to assess whether patients corresponded to monophyletic sequence groups, and the HBV genotype was determined using the genotyping tool available from the French HBV database (hbvdb.ibcp.fr). Sequences were visualized with MEGA6 and Jalview (www.jalview.org) and phylogenetic trees were edited with TreeDyn (www.treedyn.org).

Mutation rate estimation by the lethal mutation method
The number of sites that can mutate to a stop codon after a single base substitution (non-sense mutation targets, NSMTs), as well as the observed premature stop codons, were extracted from sequences. For each stop codon mutation, we identified the NSMT and calculated a correction factor C that accounts for the fact that there are three possible substitutions for each base (C = 3) and that, in some NSMTs, two of these substitutions lead to a stop codon (C = 3/2) [21]. These steps were performed using a custom R script. The preS1/preS2 and preC regions were excluded from the analysis because stop codons at these regions are not necessarily lethal and have been implicated in immune escape and pathogenesis [34,35]. For the other reading frames, we discarded sequences in which the stop codon was located in the 5% end of the protein, that lacked AUG initiation codon, with frameshift mutations, or with more than one stop codon. The mutation rate was calculated as the sum of stop codons, each multiplied by the correction factor C, divided by the total number of NSMTs.

Overall mutation spectrum
The consensus sequence of each patient was obtained and differences relative to the consensus were counted for each sequence within this patient and for each possible substitution type, using a custom R script.

Results
We analyzed 621 full-length molecular clone sequences obtained from seven studies [36][37][38][39][40][41][42], of which 451 corresponded to 22 untreated patients and 170 to 24 patients receiving long-term lamivudine/adefovir therapy (Table 1). We initially scored premature stop codons in the S, C, P, and X genes as a proxy for lethal mutations. In total, we found 23 premature stops in the S, C, P, and X genes in sequences derived from untreated patients (S1 Table).   (14), N18(44), N19 (11), N20 (26), N21 (37), N22 (18) [40] VII h Yes C T1 (9), T2 (9), T3 (8) Table 1. Below is shown the consensus sequence and the colored sequence logo summarizes variability at each position. Only variable sites are shown in the patient sequences, and stop codons are highlighted in red. As also indicated on Table 2, some stop codons appeared repeatedly in the same patient, or even in different patients. Similar alignments could be obtained for the other HBV genes by retrieving GenBank accession numbers from S1  Table 2). Based on these anomalously high counts and on previous findings suggesting that S proteins with stop codons can be maintained at high frequencies in patients [43][44][45][46], we removed this gene from subsequent analyses. Considering the counts for genes C, P, and X, the estimated HBV mutation rate was 8.7 × 10 −5 per nucleotide per cell infection cycle. This value is within the same order of magnitude as the rate reported for duck HBV and is also within the range of mutation rates reported for RNA viruses and retroviruses [18,19].
Although there were fewer available sequences for analysis in treated patients we found 21 premature stops in the C, P and X genes for this group, a significantly elevated frequency compared with untreated patients (S1 Table; Fisher exact test: p < 0.001). The estimated mutation rate for treated patients was 1.1 × 10 −3 m/n/c, i.e. 12-fold higher than in untreated patients. In the treated group, we also found that several stop codons appeared repeatedly in the same or different patients, suggesting non-lethality or the presence of mutational hotspots. Removing these cases, there were still 10 stop codons in treated patients, yielding a mutation rate estimate of 5.2 × 10 −4 m/n/c, a value still 6-fold higher than in untreated patients (Fisher exact test: p < 0.001).
These results suggest that lamivudine/adefovir treatment has a marked effect on the mutation rate of HBV in vivo. However, several possible confounders may have affected our conclusion. First, we found limited information about patient infection times, but this should not represent a major source of error because the lethal mutation rate estimation method is not sensitive to the number of viral generations elapsed. Second, some of the observed mutations might have been introduced during PCR amplification. However, for all patients except three untreated, high-fidelity polymerases with similar error rates were used (Table 1), making it unlikely that the observed differences were driven by PCR errors. Third, as mentioned above, stop-codon counts were not homogeneous among genes. However, the fact that we used always full-length genome sequences should make differences between treated and non-treated patients robust to among-gene variations. Fourth, all treated patients were infected with genotype C, whereas most sequences from untreated patients belonged to genotype B (Fig 2). Therefore, we cannot rule out the possibility that the observed differences in mutation rate were driven by a higher spontaneous mutation rate of HBV genotype C compared to genotype B. To address this issue, we downloaded 1481 molecular clone sequences of a region of the preC/C reading frames from a single study in which untreated patients infected with HBV genotype B or C were enrolled [47]. After genotyping the sequences and removing the preC region, we found four unique premature stops in the 724 sequences belonging to genotype B, versus one such stop in the 757 sequences belonging to genotype C, hence arguing against a higher spontaneous mutation rate in genotype C. Therefore, the conclusion that lamivudine/adefovir treatment increases the mutation rate of HBV appears to be robust to several sources of error and bias.
To test whether treatment modified the type of mutations produced by HBV, we scored all intra-patient single-nucleotide polymorphisms (SNPs) in the 621 full-length molecular clone sequences used above for mutation rate estimation by comparing each sequence with the consensus sequence of the patient. In total, we found 3342 intra-patient SNPs along the HBV genome, and we calculated the percentage contribution of each type of nucleotide substitution to the total number of SNPs found in each patient (mutational spectrum). The most frequent type of SNP was T!C, followed by the other base transitions A!G, C!T, and G!A, whereas base transversions were less frequent. However, the mutational spectrum differed between untreated and treated patients (Fig 3). Lamivudine/adefovir treatment was associated to a greater relative abundance of A!G SNPs (Mann-Whitney test: p < 0.001) and the complementary T!C changes (p = 0.001). All other types of nucleotide substitutions were not significantly modified in treated patients, or showed lower relative abundance.  Table 1. Sequences from untreated patients are marked in blue and those from treated patients in red. For each indicated patient, numbers show the bootstrap value of the node that delimitates the sequences from this patient. Sequences from patients N16 to N22 correspond to a study in which several family members were studied and are partially intermingled (not shown). Patients N13 and N14 are not shown because their sequences did not form welldefined monophyletic groups. Letters in the center of the tree indicate the viral genotype. doi:10.1371/journal.pone.0163363.g002

Discussion
Our spontaneous mutation rate estimate of 8.7 × 10 −5 m/n/c should be taken as an upper-limit for several reasons. First, a fraction of the observed premature stop codons may be artefacts introduced during PCR amplification, cloning, or sequencing. Second, hypothetical readthrough of some stop codons may prevent their lethality. Third, the lethal effects of some stop codons may be rescued by trans-complementation in cells co-infected with other, functional variants of the virus, leading to higher-than-expected frequencies of premature stop codons. Trans-complementation may explain the large differences among genes in the frequency of stop codons, which were >20-fold higher in the S gene than in the polymerase gene P in untreated patients (Table 2). Since reading frames are strongly overlapped in HBV, we do not interpret these variations as real changes in mutation rate along the HBV genome. We suggest that the lethality of some stop codons is more likely to be rescued by trans-complementation in the surface envelope than in the RT because HBV-infected cells produce a large excess of surface proteins, which are released to the serum in the form of empty, non-infectious sub-viral particles [48,49]. In contrast, the viral polymerase is produced in smaller amounts and should be more tightly associated with the viral genome, making trans-complementation less likely. Considering this, the mutation rate obtained for the P gene of 2.7 × 10 −5 m/n/c might be more accurate than the value obtained for other genes. Interestingly, among-gene stop-codon frequency variation was smaller in treated patients than in untreated patients (Table 2), potentially indicating a lower effect of trans-complementation. The co-infection probability increases with viral load and, hence, should tend to be higher in untreated patients than in treated patients.
Although it is commonly accepted that the main source of spontaneous mutations in HBV is the low replication fidelity of the viral polymerase, the HBV genome is subject in vivo to genome editing by cellular enzymes of the apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like 3 (A3) family, which produce in G!A substitutions in the viral genome [50,51]. Of the five premature stop codons found in X, C, and P sequences derived from untreated patients, one involved a G!A substitution (TGG!TGA), two involved C!T substitutions (unlikely to be driven by A3), and two involved base transversions (clearly not driven by A3). This suggests that A3 is not the main source of spontaneous mutations in HBV. In recent work with HIV-1, we have shown that A3 activity contributes approximately half of the spontaneous mutations found in plasma, but 98% of those found in peripheral blood mononuclear cells [33]. Here, the vast majority of sequences analyzed were obtained from plasma samples. Analysis of sequences from liver tissue may contribute to better elucidate whether A3 activity contributes significantly to the total rate of spontaneous mutation of HBV.
Lamivudine/adefovir treatment increased significantly the frequency of premature stop codons. This effect was strongest in the P gene (>15-fold) and mildest in the S gene (approximately twofold). We argue that the actual mutagenic effect of lamivudine/adefovir is best reflected by the P gene because estimation bias is probably less pronounced than in the other genes, particularly S. Our analysis of the full set of SNPs revealed a specific enrichment in A/ U!G/C substitutions in HBV sequences derived from treated patients. Recent work in mice chronically exposed to lamivudine has shown the same mutagenic pattern, as determined by sequence analysis of the highly variable D-loop region of mitochondrial DNA isolated from cortical neurons [52]. The fact that the mutational spectrum observed in treated patients was similar to that found in mouse mitochondrial DNA suggests a direct mutagenic effect of lamivudine/adefovir, as opposed to an alternative scenario in which HBV polymerase variants with modified mutational properties would have evolved in response to treatment. Since lamivudine is a cytosine analogue, it is expected to be incorporated into DNA by base-pairing with guanosine, but it is unclear how this should lead to A/U!G/C substitutions. Adefovir is an adenosine analogue and may cause A/U!G/C substitutions after being incorporated into DNA by basepairing with thymidine, followed by incorporation of cytosine by base-pairing with adefovir, or through other pathways. Further work will be required to elucidate the mechanistic basis of the observed mutagenic effects.
Our results clearly show that lamivudine/adefovir treatment increases the HBV mutation rate above 10 −4 m/n/c in each viral gene (Table 2). For a per-site mutation rate of 5.2 × 10 −4 m/ n/c, approximately 1.5 new mutations should be produced per genome after a single cell infection cycle. This suggests that, in treated patients, HBV produces mutations faster than many RNA viruses. For comparison, the spontaneous mutation rate of hepatitis C virus (HCV) in vivo estimated by the lethal mutation method is on the order of 0.3-1.2 × 10 −4 m/n/c, and ribavirin/interferon treatment increases this rate between two-and three-fold [21,32]. It is noteworthy that the effect of ribavirin/interferon treatment on the HCV mutation rate appears to be less marked than that of lamivudine/adefovir treatment for HBV. Previous work has suggested that modest increases in the mutation rate of RNA viruses can situate them close to the maximum tolerable mutation rate, or error threshold [53]. In the presence of treatments, the estimated mutation rates of HBV and HCV are similar and, therefore, we speculate that treatment may also situate HBV close to the error threshold, particularly if we consider the high degree of gene overlap shown by this virus, in which many single-point mutations modify the sequence of two proteins simultaneously. However, this extremely high mutation rate may also boost viral evolvability and promote the rapid emergence of drug-resistance and immune escape mutants. Previous work has demonstrated that, in HBV and RNA viruses, high fidelity variants can evolve in response to nucleoside analogue mutagenesis [13,[54][55][56][57][58]. If a similar process occurs in vivo for HBV, it seems not to be sufficient to compensate for the mutagenic effect of lamivudine/adefovir treatment. Future work may further elucidate how HBV evolves in response to treatment-induced mutagenesis. Ideally, this should be addressed in longitudinal studies that include samples before and after the onset of treatment, or using HBV genotypematched patients in the case of cross-section studies. Similarly, sequencing at both the acute and chronic diseases stages for patients with known infection times would allow testing whether the HBV mutation rate evolves throughout the course of infection.
Supporting Information S1 Table. List of premature stop codons indicating the patient, gene, genome position, type of stop and accession number, in Excel format. (XLSX)