Phylodynamics of HIV-1 Subtype C Epidemic in East Africa

The HIV-1 subtype C accounts for an important fraction of HIV infections in east Africa, but little is known about the genetic characteristics and evolutionary history of this epidemic. Here we reconstruct the origin and spatiotemporal dynamics of the major HIV-1 subtype C clades circulating in east Africa. A large number (n = 1,981) of subtype C pol sequences were retrieved from public databases to explore relationships between strains from the east, southern and central African regions. Maximum-likelihood phylogenetic analysis of those sequences revealed that most (>70%) strains from east Africa segregated in a single regional-specific monophyletic group, here called CEA. A second major Ethiopian subtype C lineage and a large collection of minor Kenyan and Tanzanian subtype C clades of southern African origin were also detected. A Bayesian coalescent-based method was then used to reconstruct evolutionary parameters and migration pathways of the CEA African lineage. This analysis indicates that the CEA clade most probably originated in Burundi around the early 1960s, and later spread to Ethiopia, Kenya, Tanzania and Uganda, giving rise to major country-specific monophyletic sub-clusters between the early 1970s and early 1980s. The results presented here demonstrate that a substantial proportion of subtype C infections in east Africa resulted from dissemination of a single HIV local variant, probably originated in Burundi during the 1960s. Burundi was the most important hub of dissemination of that subtype C clade in east Africa, fueling the origin of new local epidemics in Ethiopia, Kenya, Tanzania and Uganda. Subtype C lineages of southern African origin have also been introduced in east Africa, but seem to have had a much more restricted spread.


Introduction
Human immunodeficiency virus type 1 (HIV-1) sequences belonging to the pandemic group M are classified into nine subtypes (A-D, F-H, J, and K), six sub-subtypes (A1-A4, and F1-F2), and a variety of inter-subtype recombinant forms (Los Alamos HIV sequence database: http://hiv-web.lanl.gov/). Subtype C is the most prevalent variant, accounting for nearly half (48%) of all global infections [1]. This high prevalence is due to the predominance of subtype C in southern Africa, east Africa and India, with further infections in central Africa and Brazil.
Subtype C accounts for .95% of HIV infections in all southern African countries [1]. Several studies showed that subtype C sequences from neighboring southern African nations display a great degree of phylogenetic intermixing with no evidence of significant geographical clustering [2,3,4,5,6,7], indicating a largely unrestricted viral movement across the entire subcontinent. A more recent phylogenetic study revealed that after sequential pruning of ambiguously positioned taxa 10 strongly supported subtype C clusters becomes apparent in southern Africa, showing that the geographic subdivision of subtype C viruses circulating in this region is higher than expected by chance [8]. Most subtype C clusters identified, however, circulate in more than one southern African country and all four countries analyzed (Botswana, Malawi, South Africa and Zambia) comprise strains from multiple clusters. Thus, HIV epidemics in southern African countries are probably the result of the introduction and circulation of multiple subtype C strains with a variable level of local and regional dissemination.
Little is known about the genetic characteristics of HIV-1 subtype C strains circulating in east Africa. Previous studies showed that two genetically different subtype C strains designated C and C9, have been co-circulating in roughly similar prevalence and among the same risk groups and geographical areas in Ethiopia [13,15,40]. A recent study of Thomson and Fernández-García [8] revealed that the Ethiopian-C clade corresponds to one subtype C cluster also found in other east African countries including Burundi, Djibouti, Kenya, and Uganda; while the Ethiopian-C9 clade was assigned to an independent cluster associated to southern Africa. Other studies performed in Kenya showed that subtype C samples from this country are not concentrated in a single cluster, but distributed in several independent lineages associated to sequences from both east and southern Africa [34,39]. Despite these previous studies, we still have an incomplete understanding of the number, onset date, and migration pattern of the distinct HIV-1 subtype C lineages circulating in the eastern African region.
To obtain a more comprehensive picture of the spatiotemporal dynamics of the HIV-1 subtype C epidemic in east Africa, we analyzed a large number of subtype C pol sequences sampled from the east (Burundi, Ethiopia, Kenya, Tanzania and Uganda), southern (Botswana, Malawi, Mozambique, South Africa, Zambia and Zimbabwe) and central (Angola and Democratic Republic of Congo) African regions over a time period of 25 years (1986-2010).

Sequence dataset
HIV-1 subtype C pol sequences from east, southern, and central African countries, that matched the selected genomic region (nt 2253-3272 relative to HXB2 clone) were retrieved from the Los Alamos HIV Database (http://hiv.lanl.gov). Countries were grouped in geographical regions according to the classification proposed by Hemelaar et al [1]. In order to improve the accuracy of phylogenetic inference only sequences from antiretroviral therapy naïve individuals were selected. The subtype assignment of all sequences was confirmed by the REGA HIV subtyping tool v.2 [41] and by maximum likelihood (ML) phylogenetic analysis (see below) with HIV-1 subtype reference sequences. Those sequences with incorrect subtype C classification, sequences containing frame-shift mutations or deletions, multiple sequences from the same individual and those sequences from countries poorly represented (,5 sequences) were removed. This resulted in a final dataset of 1,981 HIV-1 subtype C pol sequences sampled from 12 different African countries (Table 1). Sequences were aligned using the CLUSTAL X program [42] and alignment is available from the authors upon request.

Substitution saturation and likelihood mapping analyses
Substitution saturation was evaluated by plotting the estimated number of transitions and transversions against genetic distance for each pairwise comparison in our alignment of 1,981 HIV-1 subtype C pol sequences using DAMBE program [43]. The phylogenetic signal in the pol dataset was investigated with the likelihood mapping method [44] by analyzing 10,000 random quartets. Likelihood mapping was performed with TREE-PUZ-ZLE program [45] using the online web platform Phylemon 2.0 [46].

Phylogenetic analysis
ML phylogenetic trees were inferred under the GTR+I+C 4 nucleotide substitution model, selected using the jModeltest program [47]. ML tree was reconstructed with PhyML program [48] using an online web server [49]. Heuristic tree search was performed using the SPR branch-swapping algorithm and the reliability of the obtained topology was estimated with the approximate likelihood-ratio test (aLRT) [50] based on the Shimodaira-Hasegawa-like procedure. The ML trees were visualized using the FigTree v1.3.1 program [51].

Characterization of intrasubtype C/C9 recombinant sequences
Putative intrasubtype C/C9 recombinant sequences in Ethiopia were identified by Bootscanning using Simplot version 3.5.1 [52], following the same procedure described by Pollakis et al [40]. Bootstrap values supporting branching with reference sequences were determined in Neighbor-Joining (NJ) trees constructed using the K2-P nucleotide substitution model, based on 100 resamplings, with a 300 bp sliding window moving in steps of 10 bases.

Analysis of spatiotemporal dispersion pattern
The evolutionary rate (m, units are nucleotide substitutions per site per year, subst./site/year), the age of the most recent common ancestor (T mrca , years), and the spatial dynamics of major subtype C clades from east Africa were jointly estimated using the Bayesian Markov Chain Monte Carlo (MCMC) approach implemented in the BEAST software package v1.6.2 [53,54]. Analyses were performed using the GTR+I+C 4 nucleotide substitution model, an uncorrelated Lognormal relaxed molecular clock model [55], a Bayesian Skyline coalescent tree prior [56], and a discrete phylogeographic model in which all possible reversible exchange rates between locations were equally likely [57]. Two separate MCMC chains were run for 4610 8 generations and adequate chain mixing was checked by calculating the effective sample size (ESS) after excluding an initial 10% for each run using program TRACER v1.4 [58]. MCMC runs converged to almost identical values and combined estimates showed ESS values .200. Maximum clade credibility (MCC) trees were summarized from the posterior distribution of trees with TreeAnnotator and visualized with FigTree v1.3.1. Migratory events were summarized using the cross-platform SPREAD application [59].

Phylogenetic analysis
A large dataset of HIV-1 subtype C pol sequences (n = 1,981) downloaded from the Los Alamos HIV Database (http://hiv.lanl. gov) was used to characterize the relationship between subtype C sequences sampled from east, central and southern African countries. The transition/transversion vs divergence graphics showed that both type of nucleotide substitutions increase linearly with the genetic distance, with transitions being higher than transversions ( Figure S1a), thus indicating no substitution saturation in our alignment. While, the likelihood-mapping analysis showed that most (90%) of the randomly chosen quartets from the HIV-1 subtype C alignment were equally distributed in the three corners of the likelihood map ( Figure S1b), indicating a strong tree-like phylogenetic signal in the data. Both analyses indicate that the HIV-1 subtype C pol dataset used in this study contains enough evolutionary information for reliable phylogenetic and molecular clock inferences.
The ML phylogenetic analysis revealed that most (73%) subtype C sequences from east Africa branched within a highly supported (aLRT = 0.93) monophyletic cluster, here called C EA , that contains sequences from all five east African countries analyzed ( Figure 1). Notably, the C EA clade comprises a minor proportion (9%) of the 54 sequences from central Africa, but none of the 1,576 sequences from southern Africa here included. A minor fraction (11%) of subtype C sequences from east Africa branched in a second well supported (aLRT = 0.94) monophyletic cluster that comprises sequences from Ethiopia, and corresponds to the so called Ethiopian-C9 (C9 ET ) clade ( Figure 1). The remaining subtype C east African sequences (16%) were distributed in several independent lineages of small size (n#5 sequences) that were intermixed among strains from southern African countries ( Figure 1).
The analysis of sequence distribution among clades by country of origin revealed three different patterns within east Africa represented by Burundi/Uganda, Ethiopia and Kenya/Tanzania ( Figure 2a). All or most subtype C strains circulating in Burundi and Uganda belong to the major clade C EA . Subtype C strains from Ethiopia, by contrast, were mainly distributed into clades C EA (61%) and C9 ET (37%). Finally, about 64% of subtype C sequences from Kenya and 49% from Tanzania branched within the major clade C EA , while the remaining sequences were distributed in the multiple minor clades of southern African origin. Such geographical variation in the prevalence of different subtype C clades could be also observed at a more local scale in Tanzania (Fig. 2b). In the Kagera and Mwanza regions (north), most (.70%) subtype C strains belong to the C EA clade. In the Kilimanjaro region (northeast), sequences from both the C EA and ''southern African'' clades reach a roughly similar prevalence. In the Mbeya region (southwest), only ''southern African'' clades were detected.

Migration pattern of HIV-1 C EA clade
A closer inspection of the HIV-1 C EA clade showed that sequences from Burundi occupies the most basal position in the clade ( Figure S2), thus suggesting that Burundi was the most probable epicenter of dissemination of this subtype C lineage. The migration pattern of the C EA lineage was reconstructed using a Bayesian statistical framework that allows ancestral reconstruction of the locations at the interior nodes of Bayesian tree while accommodating phylogenetic uncertainty. Sequences with no information about sampling date (n = 2), sequences with unexpectedly long branches in the phylogenetic analysis (n = 10), and Ethiopian sequences with evidence of intra-subtype recombination (n = 8, see below) were excluded from this analysis. This resulted in a final dataset of 236 sequences (Burundi = 92, Ethiopia = 47, Kenya = 24, Tanzania = 40, and Uganda = 33) sampled between 1990 and 2010.
The Bayesian MCC tree supports the hypothesis that the C EA clade originated in Burundi (PP = 1) and was later exported to the other east African countries where it further spread, establishing new local epidemics (Figures 3 and 4). Estimation of viral movement among countries, obtained by counting the state changes along the tree nodes, points to the role of Burundi as the most important hub of dissemination of this subtype C lineage in east Africa, followed by Tanzania (Table 2). Several migration events of the lineage C EA from Burundi to Ethiopia (n = 4), Kenya (n = 5), Tanzania (n = 8) and Uganda (n = 8) were detected, as well as from Tanzania to both Kenya (n = 3) and Uganda (n = 7). Importation of the C EA lineage into Burundi from other east African countries, and viral exchanges between Ethiopia, Kenya and Uganda were seldom detected in our dataset.
The Bayesian analysis also supports an important phylogeographic subdivision within the C EA lineage. Consistent with the ML topology ( Figure S2), most subtype C sequences from Ethiopia, Kenya, Tanzania and Uganda branched in countryspecific monophyletic sub-clusters that most probably (PP$0.93) had a Burundian origin (Fig. 3). The C ET1 and C ET2 lineages, that correspond to the so called Ethiopian-C clade, comprise 44% of all Ethiopian sequences here included and were almost exclusively composed by sequences from this country. The C KE and C UG lineages comprise 33% and 37% of all sequences from Kenya and Uganda, respectively, and their circulation seems to be mainly restricted to those countries. Finally, the C TZ lineage comprises 39% of all Tanzanian sequences analyzed and has also been disseminated to Kenya and Uganda. Both ML and Bayesian analyses further suggest that the C EA clade branched in two major sub-clades: one composed by sequences from Burundi and lineages C ET1 , C ET2 and C UG; the other one composed by sequences from Burundi and lineages C KE and C TZ . The statistical support of such major sub-clades in Bayesian analysis, however, was not significant (PP,0.50) and this observation should be interpreted with caution.

Time-scale of the HIV-1 C EA clade
The median estimated evolutionary rate for the pol region of the C EA clade was 1.8610 23 (95% highest posterior density [HPD]: 1.1610 23 -2.4610 23 ) subst/site/year, similar to that previously estimated for HIV-1 subtype C lineages circulating in South America [60] and southern Africa [6]. Importantly, the coefficient of rate variation was higher than zero (0.26 [95% HPD: 0.21-0.31]), thus demonstrating a significant variation of substitution rate among branches in the C EA clade and supporting the use of a relaxed molecular clock model to reconstruct the time-scale of this lineage. According to this analysis the C EA clade started to spread in Burundi at 1962 (95% HPD: 1942-1975), while major subclades C ET1 /C ET2 , C KE , C TZ and C UG began to expand in Ethiopia, Kenya, Tanzania and Uganda, respectively, by the early 1970s ( Figure 3).

Time-scale of the HIV-1 subtype C Ethiopian clades
The time-scale of the two major Ethiopian clades (C ET and C9 ET ) was also estimated by combining all sequences from this country in a single dataset and incorporating the posterior distribution of the substitution rate previously estimated for the C EA lineage as an informative prior. This analysis resulted in a Bayesian MCC tree in which clades C ET and C9 ET were poorly supported (PP,0.5) and several strains branched outside those major clades ( Figure S3). A careful exploration of Ethiopian sequences, revealed that some strains initially classified within clades C ET (n = 8) or C9 ET (n = 10) and those strains that branched outside major Ethiopian clades (n = 2) are putative C/C9 intrasubtype recombinant viruses ( Figure S3). When those viruses were excluded, the clades C ET and C9 ET segregate in two highly supported (PP.0.9) reciprocally monophyletic groups ( Figure 5). According to this new Bayesian MCC tree, the median T mrca was estimated at 1978 for clade C ET , 1981 for sub-clade C ET1 , 1984 for sub-clade C ET2 , and 1981 for clade C9 ET ( Figure 5).

Discussion
This study demonstrates a significant phylogeographic subdivision of HIV-1 subtype C strains circulating in the east respect to those circulating in the central and southern African regions, consistent with a recent study [8]. Most (73%) subtype C sequences from east Africa analyzed in this study branched within a highly supported monophyletic clade, here called C EA , that comprise 100% of subtype C sequences from Burundi, 97% from Uganda, 64% from Kenya, 61% from Ethiopia, and 49% from Tanzania. This major east African clade also comprises a minor proportion (,10%) of sequences from central Africa, but no sequence from southern Africa, thus indicating that its circulation is mainly restricted to the east African region. Of note, the genealogies previously inferred for HIV-1 subtypes A and D also support a model of limited introduction of each subtype into east Africa, followed by a subsequent local expansion [61].
Our phylogeographic study suggests that the C EA clade most probably originated in Burundi and after a period of local expansion, this viral lineage was disseminated at multiple times to Ethiopia, Kenya, Tanzania and Uganda, where it generated new local epidemics. Several introductions of the C EA lineage from Tanzania into both Kenya and Uganda were also detected, while viral exchanges between Ethiopia, Kenya and Uganda were less frequent. Five major country-specific monophyletic sub-clusters were detected within the C EA clade that comprise 44%, 33%, 37%, and 39% of all sequences from Ethiopia, Kenya, Uganda and Tanzania here included, respectively. Thus, despite frequent viral movement among east African countries, a significant proportion of subtype C infections in Ethiopia, Kenya, Tanzania and Uganda most likely resulted from the expansion of a few ancestral C EA strains.
It has been suggested that interconnectivity between population centers was a critical factor in the spread of HIV-1 subtypes A and D across Africa [61]. The restricted circulation of the C EA lineage in southern African countries is consistent with this model, considering the relative inaccessibility between the principal population centers of eastern and southern African regions. This model, however, is not consistent with the proposed role of Burundi as the main hub of dissemination of the C EA clade in the region, as this small country is poorly interconnected to other east African countries. Previous studies have also shown a strongly supported phylogenetic relationship between subtype C sequences  Table 1. The color of branches represents the geographic region from where the subtype C sequences originated, according to the map given in the figure. The boxes highlight the position of the major east African subtype C lineages. The tree was rooted using HIV-1 subtype A1 and D reference sequences (black branches). Horizontal branch lengths are drawn to scale with the bar at the bottom indicating nucleotide substitutions per site. doi:10.1371/journal.pone.0041904.g001

HIV-1 Subtype C Epidemic in East Africa
PLoS ONE | www.plosone.org from Brazil, the UK, Burundi and Kenya; thus indicating that the C EA clade has also been disseminated to South America and Europe [60,62,63]. These evidences suggest that factors other than accessibility may have shaped the dissemination of the C EA clade at both local and global scale.
Burundi has known many violent ethnic conflicts mainly since the 1960s that resulted in large migration flows. Two major civil conflicts that took place in Burundi in 1972 and 1993 generated especially large human movements with the former producing around 300,000 refugees and the latter producing about 687,000 [64]. Most refugees initially crossed the border of their country in the east, fleeing to neighboring Tanzania, followed by movement into other neighboring African countries and later to Europe and North America. It has been estimated that there are about 200,000 Burundians currently living in Tanzania, 18,000 in the Democratic Republic of the Congo, 4,000 in Uganda, 10,000 in the  European Union, and about 3,000 in the USA and Canada [64]. The molecular clock analysis clade traced the origin of the C EA lineage in Burundi to the early 1960s, while the onset date of the major sub-clades circulating in Ethiopia, Kenya, Tanzania and Uganda was estimated at around the early 1970s, coinciding with the first large Burundian migration flow. These analyses support the notion that the Burundian migration flow occurring in 1972 may have played a fundamental role in the regional and international dissemination of the C EA clade.
While subtype C epidemic in Burundi and Uganda is largely dominated by the C EA clade, a second major subtype C lineage is also circulating in Ethiopia. Our results showed that the two Ethiopian lineages previously designated C and C9 [13], resulted from independent founder strains originated in the eastern and southern African regions, respectively, and further confirmed the circulation of intra-subtype C/C9 recombinants in Ethiopia [40]. The prevalence of C/C9 recombinant viruses estimated in our dataset (20%) was equal to the percentage found in the general Ethiopian population [40]. The onset date of Ethiopian clades C and C9 was dated to between the early 1970s and the early 1980s; consistent with previous estimations [65,66,67].
A large collection of minor subtype C lineages of southern African origin were detected in Kenya and Tanzania, which together represent 36% and 51% of sequences from those countries here analyzed, respectively. These lineages seem to have a more restricted expansion than the C EA clade, although they were particularly prevalent (100%) in southwest Tanzania (Mbeya region), close to Zambia and Malawi. The co-circulation of subtype C sequences from both east and southern African origin in Tanzania is consistent with its intermediate geographical position between eastern and southern countries. It is unclear whether subtype C clades of southern African origin detected in Kenya were introduced from Tanzania and/or directly from southern Africa.
It is also unclear the relevance of these findings for HIV-1 vaccine design. Possible correlations of distinct HIV-1 subtype C clades with differential susceptibility to neutralizing antibody and/ or cellular immune responses should be explored to justify the selection of vaccines incorporating one or multiple immunogens derived from major African subtype C clades [8]. It is also uncertain whether distinct subtype C lineages may possess different biological properties that affect disease progression and viral transmission. A recent study conducted in Ethiopia showed that infection with clade C ET is associated with initially lower HIV-1 RNA plasma loads but more rapid onset of disease than infections with clade C9 ET [68]. The authors proposed that the clade C ET may be less efficiently transmitted than clade C9 ET , which is consistent with epidemiological evidence that show that the strain C9 ET has gained ground and surpassed the clade C ET over time [40,68]. New studies are necessary to determine if subtype C lineages of east African origin are less transmissible than those originated in southern Africa.
In conclusion, the results presented here point to the existence of a HIV-1 subtype C lineage characteristic of east Africa, which accounts for .70% of subtype C infections in this African region. This lineage probably emerged in Burundi in the 1960s and about 10 years later spread to Ethiopia, Kenya, Uganda and Tanzania, where it disseminated establishing new local epidemics. The subtype C epidemics in Ethiopia, Kenya and Tanzania also resulted from the introduction and dissemination of additional lineages of southern African origin. The explanation for the pattern of spread of the HIV-1 subtype C epidemic in east Africa is probably multifactorial and includes founder effects, massive migration between countries as a consequence of ethnic conflicts and geographical proximity.  years 1960, 1965, 1970, 1975, 1980, 1985, 1990 and 2000. Lines between locations represent branches in the Bayesian MCC tree along which location transition occurs. Location circle diameters are proportional to square root of the number of Bayesian MCC branches maintaining a particular location state at each time-point. The whitegreen color gradient informs the relative age of the transitions (olderrecent  Supporting Information Figure  The PP support is indicated only at key nodes. Positions of the putative interclade C/C9 recombinant sequences are marked with asterisks. Horizontal branch lengths are drawn to scale with the scale at the bottom indicating years. The tree was automatically rooted under the assumption of a relaxed molecular clock. Representative bootscanning plots of some putative C/C9 intrasubtype recombinant sequences are depicted on the right. Query sequences were compared to reference sequences of HIV-1 clades A1 (AB253429), D (AY371157), C ET (AY242589), and C9 ET (AY242581). (PPTX) Figure 5. Time-scaled Bayesian MCC tree of major Ethiopian HIV-1 subtype C lineages. MCC tree was obtained after exclusion of putative C/C9 intrasubtype recombinant sequences. Branches are colored according to the initial clade assignment of each sequence based on ML analysis: C ET (blue) and C9 ET (red). The PP support and the median age (with 95% HPD interval in parentheses) are indicated only at key nodes. Horizontal branch lengths are drawn to scale with the scale at the bottom indicating years. The tree was automatically rooted under the assumption of a relaxed molecular clock. doi:10.1371/journal.pone.0041904.g005