Tracing hepatitis B virus (HBV) genotype B5 (formerly B6) evolutionary history in the circumpolar Arctic through phylogeographic modelling

Serum samples, HBV DNA extraction and sequencing analysis

Twenty-two serum samples from various locations in Nunavut were collected from HBV/B5-infected individuals having self-identified as Inuit. All specimens were described previously in studies for which informed patient consent was obtained for HBV molecular analysis and study approval was granted from institutional ethics review boards (Health Canada/Public Health Agency of Canada research ethics board protocol number REB-2006-0048 and REB-2012-0062; University of Manitoba research ethics board approval number HS15821; Huynh et al., 2014; Larke et al., 1987; Minuk et al., 2013; Osiowy, Larke & Giles, 2011). HBV DNA was extracted from 200 µl sera by SDS-proteinase K lysis and phenol chloroform extraction methods (Osiowy, 2002) and resuspended in 30 µl sterile, nuclease-free water. Amplification of the full HBV genome was performed using a high fidelity polymerase (Roche Expand High Fidelity^Plus System, Roche Diagnostics, Laval, QC, Canada) as described previously (Osiowy et al., 2010). Samples having an HBV viral load precluding full genome amplification and sequencing were partially sequenced using primers and methods described previously (Minuk et al., 2012). Specific amplicons were gel-purified prior to cycle sequencing with an AB 3730 XL DNA Analyzer using Big Dye 3.1 terminator chemistry (Thermo Fisher Scientific, Burlington, ON, Canada). Sequences were assembled and analysed using DNA sequence analysis software (Lasergene software suite v 10.0, DNASTAR, Madison, WI, USA). Sequences were submitted to GenBank under accession numbers KP659234 –KP659255. Nine other sequences previously described from individuals residing in Nunavut were included (JN792894, JN792896, JN792897, JN792900, JN792901, DQ463795, DQ463796, DQ463799, DQ463802; Kowalec et al., 2013; Osiowy et al., 2006). Thus, a total of 31 partial and full genome HBV sequences from Nunavut were included in the analysis (Table 1).

Ten serum samples from various locations in Alaska were collected and the HBV DNA extracted and sequenced as described for Nunavut specimens. Participants were enrolled statewide by the Liver Disease and Hepatitis Program at the Alaska Native Tribal Health Consortium (Anchorage, AK). This study was approved by the Alaska Area (Indian Health Service) and the Centers for Disease Control and Prevention Institutional Review Boards (Protocol number AAIRB 2001-07-022/1996-01-001). The study and manuscript were approved by the Alaska Native Tribal Health Consortium and the Southcentral Foundation Board of Directors. All participants provided written informed consent. Sequences were submitted to GenBank under accession numbers KP659219 –KP659228. Six other sequences previously described from individuals residing in Alaska were included (AB287314 –AB297319; Sakamoto et al., 2007). Thus a total of 16 partial and full genome HBV sequences from Alaska were included in the analysis (Table 1).

Four serum specimens collected from East Greenland were extracted for HBV DNA and sequenced as described for Nunavut specimens. These specimens were collected previously for the study by Krarup et al. (2008), in which ethics approval was granted by the Commission for Scientific Research in Greenland (Approval number 505–99), and written informed consent was obtained from each participant. Sequences were submitted to GenBank under accession numbers KP659230 –KP659233. Six other sequences previously described from individuals located in West Greenland were included (AB287320 –AB287325; Sakamoto et al., 2007). Thus a total of 10 partial and full genome HBV sequences from Greenland were included in the analysis (Table 1).

Methods for specimen and data collection and analysis were carried out in accordance with the Tri-Council Policy on Ethical Conduct for Research Involving Humans (Canadian Institutes of Health Research, Natural Sciences and Engineering Research Council of Canada, and the Social Sciences and Humanities Research Council of Canada).

Phylogenetic analysis and evolutionary dynamics

All sequences listed in Table 1, including partial and full genome sequences, were aligned using ClustalX (Thompson et al., 1997), resulting in an alignment with 3220 sites. Bayesian analysis was performed using Markov Chain Monte Carlo (MCMC) methods implemented in BEAST 2 (Bouckaert et al., 2014) for the phylogeny and estimation of the HBV effective population size. In order to facilitate convergence, the MCMC chains were run sufficiently long; at 40 million generations with sampling every 10,000 steps and the first 10% of samples discarded as burn-in. During analysis, the sequences in the West Greenland clade were restricted to be monophyletic. Without monophyly constraint the analysis did not converge, and from our experience, sequences from the western circumpolar regions of Alaska, Canada and Greenland, tend to cluster into monophyletic clades (Kowalec et al., 2013; Sakamoto et al., 2007). A coalescent tree prior was chosen due to the intra-species analysis and the assumption that sampling across all clades was consistent and accurate. The XML file including all data for BEAST 2 analysis is available as Data S1.

The aligned sequence data consisting of full length and partial sequence was partitioned into eight parts consisting of sites 1–835, 836–1373, 1374–1620, 1621–1900, 1901–2307, 2308–2450, 2451–2847, and 2848–3220, based on gene boundaries in the genome. Since each of these partitions code for different genes, and some even two genes in different reading frames, these partitions will be governed by different evolutionary mechanisms. Therefore, a separate substitution model was used with each partition having its own relative substitution rate. We used the reversible-jump-based substitution model (Bouckaert, Alvarado-Mora & Pinho, 2013) with four gamma categories and invariant sites so that uncertainty in the substitution model choice is integrated out. This model jumps between the models F81 (Felsenstein, 1981), HKY85 (Hasegawa, Kishino & Yano, 1985), TAN93 (Tamura & Nei, 1993), TIM (Posada, 2003), EVS (Drummond & Bouckaert, 2015) and GTR (Tavaré, 1986) and automatically estimates the model parameters during the MCMC.

The phylogenetic tree was prepared using DensiTree implemented in BEAST 2, which has the advantage of being able to visualize the uncertainty in both node heights and topology, such that a qualitative analysis of all tree sets can be made (Heled & Bouckaert, 2013). Bayesian skyline plot analysis was used to estimate the relative HBV/B5 population size through time (Drummond et al., 2005). A coalescent tree prior is used, which is an appropriate prior for within species samples. Furthermore, it uses a non-parametric population function, so there is no commitment to a parametric population function; that is, a population history, such as an exponentially growing population or constant population, is not assumed beforehand. Default settings for the Bayesian skyline plot (five intervals, default hyper priors) were used.

Geographic analysis

A joint geographic and phylogenetic analysis of aligned partial and full genome sequences listed in Table 1 using the landscape aware geographic model (Bouckaert et al., 2012) was performed using BEAST 2. Essentially, the model uses a random walk based on successive steps in a random direction to model the spread of organisms across a landscape. The area of interest throughout the western circumpolar Arctic region that is being investigated contains many large landmasses among long stretches of coastline. The landscape aware model was configured to distinguish between sea, inland, and coastal regions and assumed fast (10-fold higher) dispersion rates along the coast compared to dispersion rates inland. Furthermore, there is a reluctance to get into water, but once in water, the rate of dispersal is very high (10-fold higher than along the coast).

A prior on the location of the root of the HBV phylogeographic summary tree was used to enforce its location on the Asian side of the Bering Strait, in order to minimize the geographical area used for the landscape aware model, and to parallel the likely HBV/B5 origin and dispersal from Asia through the Bering Strait region (Paraskevis et al., 2015); thus, the Asian HBV sequences do not impact the geography of the North American and Greenlandic HBV sequences.

During the MCMC run, locations of internal nodes are sampled. To position a branch onto the map, the sample location information is used for the start and end of a branch. To find the location of the mid-point of the branch of length t, the location in a 32 × 32 grid that maximises the probability of going from a start location to that grid point after time t∕2, multiplied by the probability of going from an end location to that grid point after time t∕2, is determined. The branch is continually split recursively until all intermediate locations are neighbouring grid points. This determines the most probable path for a branch, which can then be visualised by straight lines connecting grid points. Note that Bayesian analyses do not produce single summary trees (visualised as a blue solid line) but distributions over trees are represented. These are samples from the posterior in the form of a set of trees. Each of these trees follows a somewhat different trajectory and can be visualised individually to get an impression of the uncertainty in the distribution route as well as the area that potentially was visited. Instead of drawing each tree using an opaque line (as was done for the summary tree) a colored dot is drawn using translucency for every grid point in the 32 × 32 grid traversed by the trajectory associated with the tree. Time depth can be visualised using colour, ranging from light blue (older) to red (younger).

HBV/B5 phylogenetic analysis

A total of 72 HBV sequences were included in the phylogenetic analysis; fifty-seven HBV/B5 sequences, of which 40 were full genome sequence, and 11 HBV/B1 and 4 HBV/B2 GenBank-derived full genome sequences used as an outgroup (Table 1). HBV/B1 and B2 sequences were included in the phylogenetic analysis to delineate the ancestral foundation and putative geographic origin of HBV/B5, as HBV/B1 sequences are most phylogenetically similar to HBV/B5 sequences, while HBV/B2 sequences are more distantly related to both B1 and B5 (Sakamoto et al., 2007). Both HBV/B1 and B2 are localized to regions in Asia (Sakamoto et al., 2007) but have not been observed to circulate in Arctic regions. Location information, sample dates and the presence of the precore stop codon mutation (G1896A) for all sequences are listed in Table 1, as are the genomic regions covered for partial sequences. Due to privacy concerns given the small population sizes of the Inuit communities, detailed location information is not given for certain samples not previously described.

Figure 1 shows the relative substitution rates and preferred substitution models for the partitioned data following reversible-jump-based substitution model analysis. Since HBV has a complex genome with many overlapping open reading frames, partitioning can be expected to give a more realistic model of substitution, as each subgenomic region has a different substitution rate (Bouckaert, Alvarado-Mora & Pinho, 2013). The resulting general trends are similar to earlier findings of Bouckaert, Alvarado-Mora & Pinho (2013), such as the requirement for complex substitution models in sequence areas where genes overlap, and higher rates where there is no overlap.

The uncorrelated log-normal relaxed clock (Drummond et al., 2006) fit the data better than the strict clock, as resulting likelihoods and posterior were non-overlapping. Furthermore, the mean coefficient of correlation of 0.63 (95% highest probability density (HPD) Interval [0.52, 0.75]) suggests that the strict clock could be dismissed. The clock was calibrated using archeological evidence from human history, and since humans are the only host for HBV in Arctic regions we can assume they share the same history. Archaeological evidence relevant to a rapid Inuit (Neo-Eskimo) migration to the Eastern Arctic from regions of Alaska starting approximately 900 to 700 YBP (sigma = 50 year; Friesen & Arnold, 2008; McGhee, 2000; Raghavan et al., 2014) was used to calibrate the root of genotype B coalescence, with an external calibrator estimate of 647 to 953 YBP. This migration displaced earlier Dorset Paleo-Inuit peoples with very little evidence for genetic or cultural interaction between the two populations (Raghavan et al., 2014), thus representing a re-setting of the human population distribution in the Eastern Arctic. Based on the mean external calibrator results, the HBV evolutionary clock rate for this study was estimated at 4.1 × 10⁻⁵ substitutions per site per year (95% HPD interval [3.1 × 10⁻⁵, 5.1 × 10⁻⁵]), which is in keeping with median rates from most literature sources, and takes into account the faster rate noted for sequences from HBeAg negative individuals (Harrison et al., 2011), of which most HBV/B5-infected individuals were in the present study due to a high prevalence of the precore stop codon mutation A1896 (Table 1; Osiowy, Larke & Giles, 2011). The estimated clock rate gave a root height for genotype B5 of 902 YBP (95% HPD Interval [803–1,001]) following DensiTree visualization of the posterior distribution over the set of trees (Fig. 2). This dating is consistent with minimum tMRCA estimates in the literature for genotype B5 (Paraskevis et al., 2015) and corresponds to the external calibrator employed. The DensiTree topology was found to be well resolved with approximately two-thirds of the clades having over 95% posterior support.

Unexpectedly, the estimated dates derived for the MRCA for HBV/B5 geographic lineages appear to shift in an east to west fashion, such that the tMRCA for the Greenlandic HBV/B5 node is approximately 650 YBP, while the average tMRCA for the Canadian and Alaskan HBV/B5 nodes is approximately 590 YBP. The earlier Greenlandic B5 branch is well supported by an 88.7% probability, while there is uncertainty as to whether subsequent branching is between Nunavut (east and west) and Alaska B5 taxa (47.2% probability) or between east Nunavut and Alaska/west Nunavut B5 (27.1% probability). This demonstrates that the root of HBV/B5 is associated with Greenlandic HBV genomic sequences, thus suggesting phylogenetic evolution of the virus during westward dispersal back into Alaska (Fig. 2).

HBV/B5 population history and dispersal in the western circumpolar Arctic

Figure 3 shows the estimated effective number of HBV infections over time, where the y-axis represents the effective population size of HBV. A rapid expansion of the viral population is estimated to have occurred from approximately 400 YBP, which coincides with a rapid diversification shown in Fig. 2 within the same time frame. This coincides with archaeological evidence for a period of transformation in precontact Arctic society, during which a previously relatively homogeneous Thule society transformed into modern Inuit, starting approximately 500 YBP and continuing onwards for several centuries. This transformation is associated with population migrations, environmental shifts, and contact with Europeans (Finkelstein, Ross & Adams, 2009; Friesen, 2010; McGhee, 1994). An estimate of the current HBV/B5 population size in the order of 100,000 is consistent with current estimated host population size estimates based on population statistics (Central Intelligence Agency, 2015; Statistics Canada, 2014; United States Census Bureau, 2012), in keeping with each infected person hosting an entire HBV quasi-species population (Lauring & Andino, 2010). If the molecular clock rate is increased, this results in reduced population size estimates; thus, confidence in the rate used is provided by the current estimates being consistent with current host populations (Drummond & Bouckaert, 2015). Coincident timing of the HBV population expansion with historical events provides further confidence in the clock rate employed.

The posterior distribution over the HBV/B5 tree sets is shown in Fig. 4 projected onto the map of the western circumpolar Arctic region as a function of increasing time. An animation of the landscape aware geographic model analysis which demonstrates the entire dispersal estimation over time and from which Fig. 4 is taken is provided as a Video S1. The blue line in Fig. 4 represents the path of the maximum clade confidence tree, shown in Fig. 2, and the set of trees representing the posterior is projected onto the map as transparent coloured dots indicating uncertainty in the path of dispersal, especially over northern Nunavut island regions. Thus, the model assumes that HBV/B5 follows a “random walk” through a geographic area, but it does not infer the nature of the host population; i.e., it is neither assumed that the walk occurs through empty space nor populated space, but the virus distribution may approximate the associated (uninfected) host population. The summary tree has good posterior support for most clades (about two third of the clades in the summary tree have over 95% posterior support), except those lower in the tree in Alaska and Nunavut, and there was a clear, monophyletic separation of geographic HBV sequence clades, with 88.7% and 100% posterior support for the Alaska/Nunavut and Greenland/Nunavut/Alaska clades, respectively.

The advantage of the landscape aware model is that it allows features of the landscape, such as water, coastlines and inland areas, to be treated as having different dispersal rates. Original dispersal of HBV/B5 is estimated to have been more rapid along geographic coastal routes (Helgason et al., 2006; Maxwell, 1985; Morrison, 1999), similar to what has been described for HBV genotype C (Littlejohn, Locarnini & Yuen, 2016). The result of the landscape aware geographic model analysis suggests a rapid evolution of the HBV/B5 ancestor and spread along geographic coastal routes following its introduction into the Eastern Arctic from Alaska via Asia. As virus spread along the coast is estimated to be more rapid than inland, modelling tends to bend long branches along the coast, such as the branch from Asia into Alaska and further eastward. The modelling reconstruction shows Greenlandic viral evolution transpiring via the northern islands of Nunavut. The uncertainty of the actual dispersal into Greenland is shown by the maximum clade confidence tree vs. the posterior distribution (green dots in Fig. 4C) throughout Ellesmere Island. This uncertainty is likely due to the lack of HBV sequence coverage in the data set from the more northern regions of Nunavut. However, despite this limitation, the resulting analysis provides strong statistical support for an early dispersal route to Greenland through Ellesmere Island. Note the evident walk back towards Alaska in Figs. 4B and 4C.