Molecular Characterization of Hepatitis B Virus in People Living With HIV in Rural and Peri-Urban Communities in Botswana

1. Introduction

Hepatitis B virus (HBV) infection remains a significant global health concern particularly in Africa, which has the highest HBV prevalence globally at 5.8% [1]. In Botswana, the prevalence of hepatitis B surface antigen (HBsAg) varies by region, reaching levels as high as 22% in certain areas of the country [2]. Only three studies have reported occult hepatitis B infection (OBI) prevalence in Botswana, with a range from 6.6% to 33% in adults [2,3,4]. OBI is defined as the presence of replicative competent HBV deoxyribonucleic acid (DNA) in the blood and/or liver of individuals testing negative for HBsAg [5].

Understanding the genetic diversity of HBV within specific geographic locations is important for devising effective prevention and control strategies as HBV genotypes have clinical relevance. HBV genotypes are associated with vaccine efficacy [6], treatment response [7] , tendency of chronicity [8], HBsAg and hepatitis B e antigen (HBeAg) seroconversion [9]. Previous studies have identified HBV subgenotypes A1, D2, D3, and genotype E in Botswana [3,10,11,12,13]. Notably, the prevalence of subgenotype A1 varies among different demographic groups [3,10,11,14]. There is need for HBV surveillance and molecular characterization efforts especially in settings with high HIV prevalence and widespread antiretroviral treatment (ART) use, in which HBV drug resistance can be common [15]. Furthermore, HBV genomics are important in studying vaccine and treatment response, as well as transmission dynamics within the country.

OBI is not reported in national, regional, and global reports however it has clinical relevance. OBI was determined to be an independent risk factor for hepatocellular carcinoma (HCC) in one study [16]. Drug resistance associated mutations have been identified in participants with OBI [15]. Furthermore, infants born to mothers with positive HBsAg (HBsAg+) were diagnosed with OBI in one study [17]. One mechanism postulated to lead to the OBI phenotype is the presence of mutations that impair HBsAg detection. Some OBI associated mutations in different HBV open reading frames (ORFs), (pres1, pre2, surface, core, pre-core, X, and the polymerase domains) have been identified and studied in Botswana and South Africa [18,19]. The impact of these mutations on protein function was assessed using available online in silico tools as screening methods for potential candidates of functional in vitro studies. The Protein Variation Effect Analyzer (PROVEAN) was more accurate than other tools being studied [18,19].

In Botswana, the HBsAg positivity remains high, and the OBI prevalence is 3 to 4 times that of HBsAg positivity as prior studies have reported an adult OBI prevalence of 6.6% to 33% [2,3,4]. Therefore, the genetic diversity of HBV in both HBsAg+ and OBI participants needs to be further studied. We aimed to molecularly characterize HBV sequences from people living with HIV (PLWH) who tested positive for HBsAg and OBI, and to determine the impact of OBI-associated mutations on protein function in a cohort of PLWH in Botswana.

2. Materials and Methods

2.1. Study Population

Archived plasma samples from participants in the Botswana Combination Prevention Project (BCPP), that had previously tested positive for HBsAg (HBsAg+), and OBI (OBI+) were used [20]. Details of the BCPP study are described elsewhere [21,22]. Briefly, the BCPP study was a pair-matched cluster randomized study that enrolled 12,610 all consenting adults residing in a random sample of ~20% of households in 30 geographically dispersed villages throughout Botswana between the years 2013 and 2018. The main aim of the BCPP study was to assess if a combination of HIV prevention strategies would reduce HIV incidence at a community level compared to the standard of care. At baseline 3,596 BCPP study participants were PLWH and (83% of whom knew their HIV status) [21]. Our study was approved by the Human Research Development Committee (HRDC) at the Botswana Ministry of Health (HRDC number: 01028).

2.2. Laboratory Procedures

HBV screening has been described in our previous report [2]. Briefly, available plasma samples from 3304/3596 (91.9%) of PLWH in the BCPP cohort were screened for HBsAg and total core-antibodies (anti-HBc). HBsAg+ samples were further screened for HBeAg and immunoglobulin M core antibodies (anti-HBc IgM). HBV viral load was quantified for HBsAg+ samples with sufficient sample volume and performed in samples that tested negative for HBsAg to determine OBI prevalence using the Roche COBAS Ampliprep/Taqman Analyzer [2].

The QIAamp DNA Blood Mini kit (Qiagen, Hilden, Germany) was used to extract DNA from 200μL of HBsAg+ and OBI+ plasma samples according to the manufacturer’s protocol with a final elution volume of 30μL. HBV DNA was amplified using tiling primers adopted from Choga’s protocol [23] (Table S1). Briefly, two 10μM pools of tiling primers were prepared. A master mix for each primer pool with 5μL of DNA template was prepared, and HBV DNA was amplified using a protocol from our previous reports [15,24]. After amplification, these PCR products were combined, and library preparation followed suite [15,23,25]. The library was loaded into flow cells version R9.4.1 (Oxford Nanopore Technologies, Oxford, UK) and the GridION platform was used for sequencing.

2.3. Sequence Analyses

2.3.1. Genotypic and Mutational Analysis

Generated FASTQ files were uploaded into Genome Detective 1.132/1.133 for reference-based assembly of HBV [26] (last accessed 20 April 2023). The generated consensus sequences were downloaded and further imputed and assessed in AliView ver. 1.26 for viewing, trimming and alignment [27]. Geno2pheno (https://hbv.geno2pheno.org) (last accessed 11 December 2023), was used to assign HBV genotypes/subgenotypes. Genotypes were confirmed using phylogenetic analysis. BCPP generated sequences and reference HBV sequences from GenBank were used to construct a phylogenetic tree of the complete surface ORF using a Bayesian Markov chain Monte Carlo (MCMC) in the Bayesian Evolutionary Analysis by Sampling Trees (BEAST) v1.8.2 (BEAST Developers) program with a chain length of 100,000,000 and sampling every 10,000 generations. The analysis utilized an uncorrelated log-normal relaxed molecular clock, the Hasegawa, Kishino, and Yano (HKY) model, and the general time reversible model with gamma distributed rates of variation among sites and a proportion of invariable sites (GTR+G+I). Tracer v1.7 (BEAST Developers) was used to visualize results and confirm chain convergence. Tree Annotator v1.7.3 (BEAST Developers) was used to choose the maximum clade credibility tree after a 10% burn-in. The mutational profile of sequences generated form this study was compared to mutations in sequences generated in previous Botswana studies that were predominantly from PLWH (2009-2015) [3,10,11,28].

2.3.2. Impact of Occult-Associated Mutations on Protein Function

Sequences with a depth of >100 were used for this analysis. The Protein Variation Effect Analyzer (PROVEAN), available at http://provean.jcvi.org/index.php (accessed 24 April 2024) was used to determine the impact of occult-associated mutations on protein function Occult-associated mutations were defined as mutations identified only in OBI+ sequences and those that were overrepresented in OBI+ sequences versus HBsAg+ sequences.

3. Results

3.1. Participants Clinical Characteristics

Table 1 summarizes the clinical characteristics of participants whose plasma samples were used in the analysis. Most participants were female (66.4%) and had a median age of 43 (IQR: 36 – 49). Most participants had a low HBV viral load of <2000IU/ml (58.9%). Approximately 94.4% of participants were on antiretroviral therapy (ART) and were mostly on tenofovir disoproxil fumarate (TDF) containing regimen (61.6%). The TDF regimen also had emtricitabine (FTC) for all participants except for one who was on a dolutegravir/TDF regimen. The majority of participants had undetectable HIV-1 RNA (<40copies/mL) (88.8%). Median duration time on ART was 7 years (IQR: 4.7 – 9.9) Table 1.

3.2. Sequencing Success

Figure 1 shows the total number of sequences generated. Sequencing was attempted on 128 samples out of the 271 HBsAg+ samples and the success rate was 67.2% (86/128). Success rate for OBI+ was 29.2% (21/72). Samples that were successfully sequenced had HBV viral loads ranging from target not detected (TND) to >1.7 x 10⁸ IU/mL.

In total, 27 out of the 30 BCPP study sites contributed HBV sequences to this analysis (Figure S2). Samples from participants residing in Mmadinare, Molapowabojang and Otse (BCPP study sites) were not successfully sequenced. Most sequences were generated from participants in the Central district with 50 sequences, followed by the North-West district with 26 sequences. Kgatleng district had the least number of sequences generated (n=5) (Figure 2).

3.3. Genotypic Analysis

Figure 3 shows the clustering of the study sequences with reference sequences. All previous HBV sequences from Botswana clustered together and BCPP sequences clustered randomly (Figure 3). Overall, 93.5% (100/107) sequences were genotype/subgenotype A1, 2.8% (3/107) were D3 and 3.7% (4/107) were E. Among the HBsAg+ sequences, 93.0% (80/86) sequences were genotype/subgenotype A1, 3.5% (3/80) were D3 and 3.5% (3/80) were E. For OBI sequences, 95.2% (20/21) sequences were genotype A1 and 4.8% (1/21) were E.

3.4. Mutational Analysis

3.4.1. Escape Mutations

A total of 13 escape mutations were detected in 20% (18/90) of sequences with surface gene coverage (Table 2). _surfaceK122R had the highest frequency in 9/18 (50%) of participants with escape mutations. _surfaceT114S, _surfaceS114L, _surfaceC139R, _surfaceN146S and _surfaceC147Y are associated with impaired virion secretion. Vaccine escape mutations (_surfaceG130C, _surfaceN131T, _surfaceT121N) were identified in four participants. We also identified mutations that may impact HBsAg detection (_surfaceT118M, _surfaceC121R, _surfaceK122R, _surfaceG130C) (Table 2).

3.4.2. Comparison of Botswana reference HBV sequences and BCPP HBV sequences.

For all downstream analyses, only subgenotype A1 sequences were used as they constituted >93% of the sequences. Table 3 and Table 4 show previously generated Botswana sequences (REFERENCE-unique) and BCPP unique mutations with a prevalence of >20%. The full list of mutations is shown in Tables S3 and S4. BCPP sequences tended to have more unique mutations in all ORFs with a high prevalence (>20%) that were not observed in previous Botswana sequences. Some amino acid substitutions unique to BCPP sequences had a much higher prevalence, such as _preCV17F (53.6%), _xP33S (42.2%), and _surfaceI195M (55.7%). Among mutations that were identified in both sets of sequences, the prevalence of _pres2A7T and _pres2A11T was higher in the BCPP sequences (19.9% and 10.3% vs 5.3% for both mutations in the reference sequences). However, _pres2T38I was lower in the BCPP sequences (30.3% vs 41%). For the transcriptional transactivator protein (HBx), the prevalence of _xG22S, _xA21T and _xS46P was higher in the BCPP sequences than in the reference sequences (43.8%, 28.1%, and 67.2% vs 8.3%, 8.3% and 41.7%). In the surface protein, _surfaceN131T and _surfaceV194A had a lower prevalence in the BCPP sequences compared to the reference sequences (2.2% and 7.7% vs 37.5% and 14.8%) while the opposite was observed for _surfaceK122R (10.1% vs 2.3%).

In Table 4, we report mutations specific to the polymerase domains. It is noticeable that resistance associated mutations (_rtM204V, _rtL180M, _rtV173L) were unique to the BCPP cohort. In the same RT region, some mutations had a noticeably higher prevalence in the reference sequences as compared to the BCPP sequences. These are _rtV7A (40.0% vs 16.3%), _rtL53I (35.0% vs 11.4%), _rtH122N (35.0% vs 2.4%), _rtN332S (39.3% vs 15.9%) and _rtQ333K (44.4% vs 14.5%). In the terminal protein (TP) region, _tpH182Q had a noticeably higher prevalence in the BCPP cohort than in the reference sequences (21.7% vs 5.1%). _spY86H, _spS125N and _spS129N in the spacer domain occurred more frequently in the BCPP sequences than the reference sequences. In the RNase H domain, _RNaseHY116F was observed at a much higher prevalence in the BCPP cohort compared to the reference sequences (30.0% vs 8.3%).

3.4.3. Impact of Occult-Associated Mutations on Protein Function

For this analysis, we focused on sequences that had a depth of >100 and we identified mutations that were in sequences isolated from OBI participants only (_coreT142S, _tpE88R, _tpQ6H, _rtM250L, _preCW28L). Other mutations were overly represented in OBI+ sequences compared to HBsAg+ sequences. _surfaceV194A and _surfaceS55P appeared in 3/13 (23.1%) OBI participants each versus 1/53 (1.9%), HBsAg participant each, p-value 0.004. Using PROVEAN, three mutations were deemed deleterious, that is were deemed to affect protein function negatively: _tpQ6H, _surfaceV194A and _preCW28L.

Table 5. Impact of mutations identified only/overrepresented in participants with OBI.

ORF	Mutation	PROVEAN Prediction
Core	coreT142S	Neutral
Terminal protein	tpE88R	Neutral
	tpQ6H	Deleterious
Surface	surfaceV194A	Neutral
	surfaceS55P	Deleterious
Reverse transcriptase	rtM250L	Neutral
Precore	preCW28L	Deleterious

ORF; open reading frame, PROVEAN; Protein Variation Effect Analyzer.

Table 5. Impact of mutations identified only/overrepresented in participants with OBI.

ORF	Mutation	PROVEAN Prediction
Core	coreT142S	Neutral
Terminal protein	tpE88R	Neutral
	tpQ6H	Deleterious
Surface	surfaceV194A	Neutral
	surfaceS55P	Deleterious
Reverse transcriptase	rtM250L	Neutral
Precore	preCW28L	Deleterious

ORF; open reading frame, PROVEAN; Protein Variation Effect Analyzer.

4. Discussion

In this study, we identified HBV subgenotypes A1, D3, and E across a wide geographic area in Botswana, with subgenotype A1 representing more than 93% of all sequences. We also report the mutational profile of HBV in the Botswana population including mutations with deleterious impact on protein function in participants with OBI.

Our findings are consistent with prior HBV studies in Botswana which identified the same subgenotypes however with varying genotype prevalence [3,10,11]. We also report immune, vaccine, and diagnostic escape mutations in this population, some of which have been identified in other populations in Botswana [10,11,28]. _surfaceK122R was the most prevalent escape mutation and was identified only among participants with HBsAg positive serology. This mutation is associated with decreased HBsAg expression and HBsAg detection failure [33,34,35], however it was not detected among the OBI samples. We did not perform quantitative HBsAg ELISA which could have revealed HBsAg levels in samples with these mutations compared to those without. We also identified known mutations (_surfaceT118M, _surfaceN146S), and uncharacterized mutations at positions associated with immune escape (_surfaceC121R, _surfaceQ129C, _surfaceG130C,) in OBI participants although these mutations are not unique to individuals with OBI in other studies [30,41,45,46]. All these mutations were identified in the major hydrophilic region (MHR) (position 99 to 169) of the HBsAg with the majority of these being found in the “α” determinant of the MHR (position 124 to 147) which is a major cluster of antigenic epitopes [47]. We also identified vaccine escape mutations (VEMs) in 16.7% of participants with escape mutations, which is a cause for concern as these may counter-act vaccination efforts in the country. This was at position 131 of the surface region also in the “α” determinant of the MHR known for mutations that allow the virus to evade vaccine induced immune response [48]. We have identified a potential vaccine escape mutation, G130C which was reported as a novel mutation in 2017 [49]. Mutations at position 130 are reported to be vaccine escape mutations and have been isolated in vaccinated individuals [38,39].

There was a change in mutation patterns between sequences previously generated in Botswana (2009-2015) [3,10,11,28] and sequences we generated in the BCPP study (2013–2018). For example, the RT region of BCPP sequences harbors more drug resistance associated mutations than the reference sequences. Over 90% of PLWH among the BCPP participants were ART-experienced while in the previous studies, participants were mostly ART-naïve [3,10,11,28]. Due to the overlapping pattern of the HBV genome, some of the mutations observed in the RT region of the polymerase affected the surface region. Position _surfaceI195 corresponds to the _rtM204 [50], therefore its prevalence is higher in BCPP sequences compared to previous Botswana sequences. Furthermore, _surfaceE164D is known to alter HBsAg antigenicity and tends to occur with _surfaceI195M, as observed in our study and a previous study [51].

Most of the mutations identified in our study are uncharacterized. Other mutational variations of interest between BCPP and previously generated Botswana sequence are the _xP33S in the X region. This mutation was only observed in the BCPP sequences and had not been previously identified in Botswana. It is a B-cell epitope mutation that has been shown to result in increased endoplasmic reticulum (ER) stress [52] and reduced protein stability in combination with other mutations [53]. There were some mutations that were common to both sequence datasets; however, they were more prevalent in the BCPP sequences. For example, the _xG22S, which is reported to be an HCC-related HBx mutation [54]. BCPP sequences also had the _xT36A, a functionally characterized HCC-associated mutation at >20% prevalence. This mutation is reported to enhance viral genome integration into the host cell resulting in insertion mutations and a 3’-terminal truncation of HBx [55,56]. While previously generated Botswana sequences and the BCPP sequences were generated from different parts of the country, we cannot rule out the possibility that more people are progressing to chronicity and HCC in the population.

To further elucidate on OBI-associated mutations, we used a freely available online tool, PROVEAN, previously shown to be more accurate in predicting the impact of functionally characterized HBV mutations on protein function than other prediction tools [18,19]. This analysis also allows for the selection of mutations that could be candidates for further in vitro studies. Three OBI associated mutations were deemed deleterious: _tpQ6H, _surfaceV194A and _preCW28L. The _tpQ6H mutation has not been characterized however it falls within the N-terminal helices of the TP. Mutations and deletions in this subdomain were shown to impact RNA packaging, DNA synthesis and protein priming [57]. _preCW28L has not been characterized however a stop codon at this position has been identified in a participant with OBI [58] and in patients with HCC in a much older study [59]. Other surface mutations (_surfaceS55P and _surfaceV194A) were overrepresented in OBI+ sequences compared to HBsAg+ sequences, and other studies also show that these mutations are not unique to OBI+ sequences [30,60]. It is worth exploring the HBsAg levels in HBsAg+ and OBI+ samples with these mutations. _surfaceS55P has not been functionally characterized, while _surfaceV194A is a well know mutation associated with decreased extracellular HBsAg levels [61].

Our analysis focused on the differences between sequences from HBsAg+ and OBI+ infection, but it could be that the host factors have a key role to play on the OBI phenotype. Host factors such as immune response, human leukocyte antigen classA2 and interleukin-10 have been associated with OBI persistence [31,62,63]. Some studies attribute this phenotype to viral factors such as epigenetic control mechanisms [64,65] and mutations [18,19,66] that lead to reduced HBsAg expression and DNA replication..

A strength of the study is that we used samples from the BCPP study which recruited participants from different communities in Botswana. However, we analyzed sequences from only 27 out of the original 30 BCPP study sites due to insufficient volumes and the low sequencing success rate. Sequencing success was generally low, especially for OBI+ samples, however, this study still provides the largest number of HBV sequences in Botswana. We generated HBV sequences from participants with a target not detectable viral load result, a group usually excluded during sequencing. Escape mutations, _surfaceT114S, _surfaceN131T and _surfaceK122R were identified in these participants which shows that there may be an underrepresentation of mutations in different reports that may not sequence participants with these viral load results. There is need however, to adopt a nucleic acid testing assay with a lower limit of detection than the one used in our current study. A limitation of the study is that participants were all PLWH predominantly on ART which limits generalizability to the general Botswana population. We make comparisons between previously generated Botswana sequences and the BCPP sequences, which comes with a limitation as these datasets are not from the same communities. Therefore, we cannot rule out the impact of host-genetic factors which could limit the conclusions in the differences we note. However, the mutational profiles of these sequences are fully documented in this study.

5. Conclusions

We have molecularly characterized HBV from PLWH in Botswana and identified that subgenotype A1 is predominant countrywide. We have also reported vaccine escape mutations which shows the importance of periodic monitoring of circulating HBV strains in the population. It is essential to generate HBV sequencing data to monitor the evolution of HBV and the emergence of mutations that could evade immunity and vaccines, potentially affecting HBV prevention and management strategies.