Origin and Spread of HIV-1 Subtype B Among Heterosexual Individuals in Bulgaria

Abstract

Human immunodeficiency virus (HIV) was originally introduced in Bulgaria through heterosexual transmission (HET) and later transferred to other vulnerable groups along with numerous more recent introductions from outside Bulgaria. To define the diversity, origins, and dynamics of the HIV-1 subtypes prevalent in HET population in Bulgaria, we applied phylogenetic and phylodynamic analyses using polymerase (pol) sequences from HET individuals to infer the spatiotemporal evolutionary history of the HIV-1 epidemic in this population in Bulgaria. High genetic diversity was found, including 13 different HIV-1 subtypes: 45.7% subtype B, 19.9% CRF01_AE, 7.5% CRF02_AG, 7.5% sub-subtypes A1 and A6, 7.1% subtype C, 5.3% subtype F1, 4.0% URFs, 1.2% CRF05_DF, 0.6% subtype G, 0.3% CRF04_cpx, 0.3% CRF29_BF, 0.3% CRF14_BG, and 0.3% subtype H. The estimated root of the subtype B in the phylogenetic tree dated back to the year 1980 largely due to multiple introductions of subtype B from outside the country. Several significant clades have been identified highlighting six different main epidemic entrances of subtype B dating from 1989 to 2007. The Bayesian skyline plot showed two different exponential growth periods starting in the 1980s to 1990 followed by a constant phase up to about 2008, with another exponential growth period from 2008 to the year 2012. The migration analysis identified dynamic pattern of gene flow and demonstrated that many HET probably acquired the infection abroad (14.6%), while only (6.6%) of non-HET were infected outside country. The phylogenetic analysis showed an intermixing between sequences from Bulgarians with sequences from other countries, suggesting different HIV introduction in this country followed by the internal spread through local transmission networks.

Introduction

The global human immunodeficiency virus (HIV)-1 pandemic was reported first in men who have sex with men (MSM) in 1981 as the syndrome now known as AIDS.¹ Later it was found that MSM were one of the main vulnerable populations in which HIV-1 was quickly disseminated across North America and West Europe.² In Bulgaria, HIV-1 was found for the first time in 1986 in a hemophiliac individual, who repeatedly had received coagulation factors. Furthermore, in 1987–1988 a large-scale national campaign for HIV testing was conducted, and a number of people from harbor cities, large collectives, and groups vulnerable to blood borne infections were tested, resulting in diagnosis of HIV-1 infection in several dozen individuals.

Over the next 20 years 575 individuals infected with HIV-1 were identified, and in 2005, 89.0% of them reported that they were heterosexual transmission (HET). In the following years, the rate of newly registered persons who inject drugs (PWID) and MSM increased, and in 2012 the total cumulative HIV diagnoses reached 1,606, and only 56.7% of infected reported that they acquired the infection through heterosexual contact.^3
–5 Our epidemiological data indicated that in Bulgaria, HIV was originally introduced and disseminated through HET in the general population and later transferred to other vulnerable groups, while in North America and Western Europe since the beginning the epidemic was driven primarily by MSM.²

Our previous studies have shown that various HIV-1 subtypes and recombinant forms have been introduced and disseminated unequally among individuals from different transmission groups in the country, including HET, MSM, and PWID.^3

–8

Given the observed characteristics in the history of HIV-1 epidemic, defining the diversity, origin, and dynamics of the dominating HIV-1 subtypes in the HET population in Bulgaria is of great epidemiological importance. Nonetheless, little is known about how and where the main HIV-1 strains originated in HET, or if this group will be a generator for dissemination of the epidemic in other risk groups as has been reported in Eastern Europe.⁹

In this study, we performed molecular epidemiological, phylodynamics, and migration analyses using polymerase (pol) sequences from HET individuals to infer the spatiotemporal evolutionary history of subtype B, which is the most prevalent subtype of the HIV-1 epidemic in this population in Bulgaria.

Materials and Methods

Ethics statement

This study was approved by the Ethics Committee at the National Centre of Infectious and Parasitic Diseases, Sofia, Bulgaria (NCIPD IRB IORG0006384).

Study population and specimen preparation

In this nationwide study, 322 samples from HIV-1 infected heterosexual patients were collected between 2002 and 2012 at the National Reference Laboratory of HIV and/or in the five clinics responsible for the management of patients with HIV in Bulgaria, including the cities of Sofia, Plovdiv, Varna, Pleven, and Stara Zagora. The sequences were linked to demographic and clinical data through an anonymous numerical code according to the standards of the local ethics committee, and information on geographic region and date of sampling was retained.

Sequence analysis and data set

Viral RNA was extracted from plasma samples using the QIAamp^® UltraSens™ Virus Kit 50 (Cat. No. 53704; QIAGEN). The HIV-1 pol gene was sequenced using the ViroSeq HIV-1 Genotyping Test (Abbott) and/or TruGene DNA Sequencing System (Siemens Healthcare) and either the Applied Biosystems 3130xl genetic analyzer or an OpenGene DNA sequencing system following the manufacturer's protocol.⁵ HIV-1 subtypes were determined using the automated tool COMET¹⁰ and phylogenetic analysis as previously described.³ All 23-drug resistance mutation codons were manually removed from the alignment to exclude the possibility of convergent evolution at these sites.

Two different data sets were built for the expanded phylogenetic analysis. The first data set contained 125 HIV-1 B subtype Bulgarian pol gene sequences randomly selected and representative, plus 250 HIV-1 B subtype pol gene sequences from different countries (Armenia AM, Cyprus CY, Greece GR, Hungary HU, Romania RO, Russia NU, Serbia RS, Turkey TR, Ukraine UA, Montenegro ME, Italy IT, Albania AL, Austria AT) downloaded from the HIV Los Alamos database (www.hiv.lanl.gov/content/index) to investigate whether the Bulgarian clusters were monophyletic.

The second data set contained only the 125 HIV-1 B Bulgarian isolates. It was used to estimate the mean evolutionary rate, to investigate the temporal origin of the clusters by Bayesian molecular clock analysis, and to perform the population dynamics. Moreover, this data set was also used to test the hypothesis of restricted gene flow among distinct HIV-1 geographic regions in Bulgaria. The demographic characteristics for Bulgarian sequences are presented in Table 1a and b.

Table 1.

(a) and (b): Demographic and Epidemiological Characteristics of the Bulgarian Study Population

(a)
Patient's characteristic	Total HET with HIV-1, n (%)	Total non-HET with HIV-1, n (%)	Total individuals with HIV-1 in Bulgaria
Total	911 (56.73)	695 (43.28)	1606
Age (years)
<19	32 (3.5)	76 (10.9)	108 (6.73%)
20–44	709 (77.8)	589 (84.7)	1298 (80.80%)
>45	170 (18.7)	30 (4.3)	200 (12.45%)
Gender
Male	572 (62.8)	635 (91.4)	1207 (75.16%)
Female	339 (37.2)	60 (8.6)	399 (24.84%)
Likely country of infection
Bulgaria	778 (85.4)	649 (93.4)	1353 (84.25%)
Other country	133 (14.6)	46 (6.6)	253 (15.75%)
Other epidemiological data
Sex worker	19 (2.1)	35 (5.0)	54 (3.36%)
Patient with sexual transmitted infection	30 (3.3)	11 (1.6)	41 (2.55%)
Pregnant women	46 (5.0)	4 (0.6)	50 (3.11%)
Blood donors	93 (10.2)	9 (1.3)	102 (6.35%)
Diagnosis period
1986–1995	110 (12.1)	34 (4.9)	144 (8.97%)
1996–2004	313 (34.4)	35 (5.0)	348 (21.67%)
2005–2012	488 (53.6)	626 (90.1)	1114 (69.36%)

(b)
Patient's characteristic	HET with HIV-1 sequences available, n (%)
Total	322 (35.3) of HET
Age (years)
<19	13 (4.0)
20–44	254 (78.9)
>45	55 (17.1)
Gender
Male	191 (59.3)
Female	131 (40.7)
Likely country of infection
Bulgaria	290 (90.1)
Other country	32 (9.9)
Other epidemiological data
Sex worker	7 (2.2)
Patient with sexual transmitted infection	11 (3.4)
Pregnant women	24 (7.5)
Blood donors	32 (9.9)
Diagnosis period
1986–1995	18 (5.6)
1996–2004	90 (28.0)
2005–2012	214 (66.5)

HET, heterosexual transmission; HIV, human immunodeficiency virus.

Likelihood mapping analysis

The phylogenetic signal of each sequence data set was investigated by means of the likelihood mapping analysis of 10,000 random quartets, generated using TreePuzzle.¹¹ For a quartet, just three unrooted tree topologies are possible. The likelihood of each topology was estimated with the maximum likelihood (ML) method, and the three likelihoods are reported as a dot in an equilateral triangle (the likelihood map). Three main areas in the map can be distinguished: the three corners representing fully resolved tree topologies, that is, the presence of tree-like phylogenetic signal in the data; the center, which represents star-like phylogeny, and the three areas on the sides indicating network-like phylogeny, that is, presence of recombination or conflicting phylogenetic signals. When using this strategy, if more than 33% of the dots fall into the center of the triangle, the data are considered unreliable for the purposes of phylogenetic inference.

Phylogenetic analysis

For each data set, sequence alignments were obtained with ClustalW (BioEdit software) followed by manual editing, as already described.¹² The evolutionary model was chosen for each data set as the best-fitting nucleotide substitution model in accordance with the results of the hierarchical likelihood ratio test implemented in Modeltest software version 3.7, as already described.¹³

ML phylogenetic tree of the first data set was built with FastTree using GTR+I+G.¹⁴ The statistical robustness and reliability of the branching order within the phylogenetic trees was confirmed by the local support values (Posterior Probability >0.90).

Evolutionary rate estimates, time-scaled phylogeny reconstruction, and demographic history

The mean evolutionary rate was estimated on the second data set (HIV-1 B subtype pol gene isolates from Bulgaria) using a Bayesian–Markov–Chain–Monte–Carlo (MCMC) method implemented in BEAST package v. 1.8.0 implementing the GTR+Invariant+Gamma model using both a strict and an uncorrelated log-normal relaxed clock model.¹⁵ As coalescent priors, four parametric demographic models of population growth (constant size, exponential, logistic growth, expansion) and a Bayesian skyline plot (BSP, a nonparametric piecewise constant model) were compared.

Chains were conducted for at least 150 × 10⁶ generations and sampled every 15,000 steps. Convergence was assessed on the basis of the effective sampling size (ESS) using Tracer software v. 1.5 (http://tree.bio.ed.ac.uk/software/tracer/). Only parameter estimates with ESSs of >200 were accepted. Uncertainty in the estimates was indicated by 95% highest posterior density (95% HPD) intervals, and the best fitting models were selected by a Bayes factor (BF, using marginal likelihoods) implemented in Beast.¹⁶ In accordance with Kass and Raftery,¹⁷ the strength of the evidence against H0 was evaluated as follows: 2lnBF <2 no evidence; 2–6 weak evidence; 6–10 strong evidence, and >10 very strong evidence. A negative 2lnBF indicates evidence in favor of H0. Only values of >6 were considered significant.

To reconstruct the time-scaled phylogeny on the second data set with a Bayesian MCMC method, the molecular clock, the demographic model, and the mean evolutionary rate, previously estimated, were used. The MCMC chains were run for at least 150 million generations and sampled every 15,000 steps. Convergence was assessed on the basis of the ESS, and only parameter estimates with ESSs of >200 were accepted. Statistical support for specific clades was obtained by calculating the posterior probability of each monophyletic clade. Uncertainty in the estimates was indicated by 95% HPD intervals. The demographic history was analyzed on the second dataset by performing the BSP.

Viral gene flow analysis

The MacClade version 4 program (Sinauer Associates, Sunderland, MA) was used to test viral gene out/in flow among HIV-1 B subtype infected subjects (second data set) using a modified version of the Slatkin and Maddison test as already described.¹⁸ A one-character data matrix was obtained from the original data set by assigning to each taxon (viral sequence) in the tree a one-letter code indicating different places of Bulgaria. The putative origin of each ancestral sequence (i.e., internal node) in the tree was inferred by finding the most parsimonious reconstruction (MPR) of the ancestral character. The final tree length, that is, the number of observed viral gene flow events in the genealogy, can easily be computed and compared to the tree-length distribution of 10,000 trees obtained by random joining-splitting (null distribution). Observed genealogies significantly shorter than random trees indicate the presence of subdivided populations with restricted gene flow.

The viral gene flow among different places (character states) was traced with the State changes and stasis tool (MacClade software), which counts the number of changes in a tree for each pairwise character state. When multiple MPRs were present (as in our data sets), the algorithm calculated the average migration count over all possible MPRs for each pair.

Statistical analysis

Fisher's exact test was used to compare the demographic and epidemiologic characteristics of HIV-1-infected HET and other risk groups.¹⁹ Two-tailed tests were used to determine the significance of the comparisons using the software GraphPad Prism v5.0.²⁰

Results

Study population demographics and HIV-1 prevalence by risk group

The demographic and epidemiological characteristics of HET individuals with HIV-1 were compared to those from the other HIV-1-positive populations. As shown in Table 1a, 62.8% of HIV-1-infected HET were male, with almost equal prevalence of male individuals among HETs compared to other risk groups, including MSM and PWID (62.8% vs. 91.4%, p < .01). Individual ages of HETs varied between 17 and 67 years, with a mean of 33.4 years. The proportion of youth individuals (below 19 years) was significantly lower among HETs compared to other transmission groups (3.5% vs. 10.9%, p < .01) and contrary to the proportion of the adult, group (over 45 years) was significantly higher in HETs compared to MSM and PWID (18.7% vs. 4.3%, p < .01).

The vast majority of the patients were infected in Bulgaria (Table 1a); however, more HET individuals indicated infection abroad compared to other transmitted groups, including MSM and PWID (14.6% vs. 6.6% p < .01). The proportion of HET with HIV-1 sequences available was 35.3% (n = 322) (Table 1b). Furthermore, HETs were more evenly distributed in the country, including the capital city of Sofia, three other big cities, as well as other regions, compared to PWID that were reported in our previous study to be concentrated mostly in Sofia and Plovdiv and affected by specific CRF (Table 2).

Table 2.

HIV-1 Subtype Distribution Among Tested Population

Subtype	Number (%)	Male, n (%)	Female, n (%)	Migrants, n (%)	Sofia, n (%)	Plovdiv+ Varna+ Burgas, n (%)	Other regions, n (%)
B	147 (45.7)	105 (55.0)	42 (32.1)	18 (38.3)	72 (55.0)	36 (43.4)	39 (35.8)
01_AE	64 (19.9)	29 (15.2)	35 (26.7)	0	29 (22.1)	9 (10.8)	27 (24.8)
02_AG	24 (7.5)	9 (4.7)	15 (11.5)	4 (8.5)	9 (6.9)	10 (12)	5 (4.6)
A1 and A6	24 (7.5)	12 (6.3)	12 (9.2)	6 (12.8)	6 (4.6)	6 (7.2)	12 (11.0)
C	23 (7.1)	14 (7.3)	9 (6.9)	6 (12.8)	4 (3.1)	6 (7.2)	13 (11.9)
F1	17 (5.3)	10 (5.2)	7 (5.3)	4 (8.5)	5 (3.8)	6 (7.2)	6 (5.5)
Unclassified	13 (4.0)	9 (4.7)	4 (3.1)	5 (10.6)	4 (3.1)	6 (7.2)	3 (2.8)
05_DF	4 (1.2)	0	4 (3.1)	0	0	3 (3.6)	1 (0.9)
G	2 (0.6)	1 (0.5)	1 (0.8)	2 (4.3)	2 (1.5)	0	0
04_cpx	1 (0.3)	0	1 (0.8)	0	0	0	1 (0.9)
29_BF	1 (0.3)	0	1 (0.8)	1 (2.1)	0	0	1 (0.9)
14_BG	1 (0.3)	1 (0.5)	0	1 (2.1)	0	1 (1.2)	0
H	1 (0.3)	1 (0.5)	0	0	0	0	1 (0.9)
Total	322	191	131	47	131	83	109

HIV-1 diversity and distribution of subtypes in HET

The phylogenetic analysis revealed a broad distribution of various subtypes in the subgroup of HET (n = 322), including significant numbers of CRFs, unclassified URFs, and other subtypes of HIV-1. Of the 322 studied sequences, 147 (45.7%) were classified as subtype B, 64 (19.9%) CRF01_AE, 24 (7.5%) CRF02_AG, 24 (7.5%) sub-subtypes A1 and A6, 23 (7.1%) subtype C, 17 (5.3%) subtype F1, 13 (4.0%) URFs, 4 (1.2%) CRF05_DF, 2 (0.6%) subtype G, 1 (0.3%) CRF04_cpx, 1 (0.3%) CRF29_BF, 1 (0.3%) CRF14_BG, and 1 (0.3%) subtype H (Table 2). Analysis of the gender distribution of viral genotypes revealed a significantly higher prevalence of subtype B in male 105 (71.4%) compared to female individuals with only 42 (28.6%) while the other subtypes and CRFs were distributed more evenly between both genders. Geographic distribution showed that most of the infections with the main subtype B were found in individuals from the capital city of Sofia 72 (49%), and the rest were distributed across the entire country. CRF01_AE was predominant in Sofia with 26 (45.3%) followed by other regions, while CRF02_AG was most prevalent in Plovdiv with 10 (41.7%).

Likelihood mapping and phylogenetic analysis

The phylogenetic signal of each data set was investigated by means of likelihood mapping. The percentage of dots falling in the central area of the triangles was 24.2% and 18.3% for the first and second data set, respectively (data not shown). ML phylogenetic tree of the first data set is reported in Figure 1 and Supplementary Table 1; Supplementary Data are available online at www.liebertpub.com/aid. The tree shows that, although sequences from other countries are sometimes intermixed with the Bulgarian isolates, several statistically significant clusters, including only sequences from Bulgaria, can also be detected. Bulgarian sequences that cluster together were isolated from individuals originated from different cities of the country. Five main statistically supported clusters (cluster 6, 8, 10, 12, and 17), including only isolates from Bulgaria, are highlighted in Figure 1. Cluster 6 included infected individuals with a median age of 25 years, and 66.7% were female (F). All the individuals were residents in Sofia; moreover, epidemiological information suggested Bulgaria as probable country of infection, this could also be supported by the monophyletic characteristic of this cluster. Cluster 8 was composed of infected individuals with a median age of 29 years, and 10% were male (M). The 80% of the individuals included in this cluster were residents in Sofia and 20% at St. Zagora. Cluster 10 included infected individuals with a median age of 26 years, and 82.4% were female. This cluster included 34 individuals with a heterogeneous situation regarding residence (79.5% Sofia, 2.9% Lom, 6.0% Pazardzhik, 2.9% Razgrad, 2.9% Shumen, 2.9% St. Zagora, 2.9% Vidin). The 14.7% of the infected individuals in this cluster were represented by Gypsy population, 23.5% blood donor, 5.9% had also other sexual transmitted infections, and 5.9% as prisoner. Moreover, epidemiological information available reported contact between two patients (953 and 524), which clustered together in the phylogenetic analysis. Cluster 12 which included infected individuals with a median age of 37 years was characterized by 62.5% of female. The 87.5% of the individuals were residents in Varna and the 12.5% at Ruse. Epidemiological data reported contact between two patients (351 and 355) which clustered together and were linked in the phylogenetic tree. Finally, cluster 17 included 100% F with median age 32. A heterogeneous situation regarding residence was also showed for this cluster (44.5% Sofia, 11.1% Burgas, 22.2% Plovdiv, 11.1% Dobrich, 11.1% Varna).

FIG. 1.

Maximum likelihood phylogenetic analysis of the first data set. The tree was rooted by the midpoint rooting. Branch lengths were estimated with the best fitting nucleotide substitution model according to a hierarchical likelihood ratio test and were drawn to scale with the bar at the bottom indicating 0.02 nucleotide substitutions per site. One asterisk (*) along the branches represents significant statistical support for the clade subtending that branch (local support values >90%). The main cluster has been collapsed and highlighted. The sequences from different countries were represented by an acronym (Romania = RO; Russia = RUS; Italy = IT; Turkey = TR; Cyprus = CY; Greece = GR; Serbia = RS; Montenegro = ME; Ukraine = UA; Albania = AL; Hungary = HU; Austria = AT; Armenia = AM; Bulgaria = BG). The gender and age were shown close to the highlighted cluster.

Evolutionary rate estimates, time-scaled phylogeny reconstruction, and demographic history

The mean evolutionary rate was estimated in 125 isolates of the most prevalent HIV-1 B subtype in Bulgaria (second data set). Comparison of the marginal likelihoods of the trees showed that the relaxed clock model fitted the data better than the strict clock (2lnBF = 70.8).

Under the relaxed clock, the BF analysis showed that the BSP was better than the other models (2lnBF >89). The estimated mean value of the HIV-1 B subtype pol gene evolutionary rate of the second data set was 1.7 × 10⁻³ substitution/site/year (95% HPD: 1.026 × 10⁻³–2.37 × 10⁻³).

Figure 2 shows the Bayesian maximum clade credibility tree of the second data set and the time of the most recent common ancestor (tMRCA) estimation. The root of the tree dated back to the year 1980 (95% HPD: 1975–1986). Several significant clades have been identified. Specifically six main clades (I, II, III, IV, V, VI) and two subclades (VI1 and VI2) were highlighted. Clade I included four sequences (three collected from female and one from male) and dated back to the year 1999 (95% HPD: 1989–2001). Clade II, which included nine sequences (from female individuals), dated back to the year 2004 (95% HPD: 2002–2009). Clade III was composed of nine sequences, five of which were collected from female individuals and four from male, dated back to the year 1989 (95% HPD: 1986–1998). Clade IV which included four sequences (three from female and one from male) dated back to the year 1998 (95% HPD: 1996–2006). Clade V was composed of five sequences collected from male individuals and dated back to 2007 (95% HPD 2004–2010). Clade VI (which included 35 sequences) dated back to the year 1989 (95% HPD: 1986–1996). Inside clade VI, two statistically supported subclades (VI1 and VI2) were identified. Subclade VI1 dated back to the year 1997 (95% HPD: 1992–2000); subclade VI2 which included six sequences (five isolated from female and one from male individuals) dated back to the year 2003 (95% HPD: 2001–2007).

FIG. 2.

Bayesian phylogenetic tree of the second data set. The time of the most recent common ancestor, with the credibility interval based on 95% HPD interval, was reported in years. One asterisk (*) along the branches represents significant statistical support for the clade subtending that branch (posterior probability >0.98). Scale years are reported at the bottom of the figure. HPD, highest posterior density.

The BSP (Fig. 3) showed an exponential phase from 1980s to 1990 followed by a constant phase up to about 2008; another exponential growth from 2008 to the year 2012 can be identified followed by a constant phase.

FIG. 3.

Bayesian skyline plot of the second data set. The effective number of infections is reported on the Y-axis. Time is reported on the X-axis. The gray area corresponds to the credibility interval based on 95% HPD interval.

Viral gene flow analysis

HIV-1 B subtype Bulgarian sequences were assigned to four geographical groups (Northeast; Northwest; Southeast; Southwest) corresponding to the main cities of the country, respectively: Varna, Montana, Burgas, Sofia, and Plovdiv, and gene flow among the groups was estimated (Fig. 4). The null hypothesis of panmixia (i.e., no population subdivision or complete intermixing of sequences from different geographic areas) was rejected by the randomization test (p < .0001). Most of the observed viral gene flow events (31.6%) occurred from Southwest to Northeast (Fig. 4). A 10.5% of gene flow was also observed from Southeast to Southwest (Fig. 4).

FIG. 4.

Migration pattern of HIV-1 B subtype circulation in Bulgaria. The bubblegram shows the frequency of gene flow (migrations) in Bulgaria to/from different geographic areas. The surface of each circle is proportional to the percentage of observed migrations in the maximum likelihood genealogy. Migrations were inferred on the second data set with a modified version of the Slatkin and Maddison algorithm (1): Northeast; (2): Northwest; (3): Southeast; (4): Southwest. HIV, human immunodeficiency virus.

Discussion

In this nationwide study, we analyzed the origin, spread, and characteristics of the dominating strain of the HIV-1 subepidemic among HET populations in Bulgaria. As previously reported, the early HIV-1 epidemic in Bulgaria was driven mostly by introduction and dissemination of the virus by HET with PWID and MSM composing only a small fraction of infections.⁶

We found that more than half of the HET individuals diagnosed with HIV-1 were men and yet the dissemination between genders was more evenly distributed compared to non-HET transmission groups MSM and PWID, which is in line with our previous study.⁵ However, in this group men outnumbered women, which may indicate that some male individuals may have sexual contacts with men, or both women and men, but they have not indicated it in the self-reporting questionnaire during diagnostic procedures.

The youngest individuals with HIV-1 less than 19 years were significantly less in HET populations than in other vulnerable groups. On the contrary, the eldest individuals were more often found in HET compared to MSM and PWID. These data affirm that teenagers are at increased risk of acquiring HIV-1 infection if they participate in a vulnerable group.

Our epidemiological data indicated that significantly more HET individuals reported that they acquired the infection abroad, while fewer non-HET were infected outside country. This epidemiological characteristic indicates that there is a stream of viral strains in HET introduced from different countries of the world leading to a variety of HIV-1 subtypes in the general population. Furthermore, these strains could spill over into other transmission groups, including MSM and PWID, causing rapidly growing outbreaks in most of these risk populations that were already reported in our previous study.⁵

The assessment of viral diversity, epidemiological history, and migration pathway are key factors in monitoring of the HIV-1 epidemic.

Our previous studies^3

–6 have focused to analyze the key characteristics of the HIV-1 epidemic in Bulgaria highlighting high genetic diversity, including subtype B representing the dominant clade and the presence of a number of CRFs, URFs, and rare HIV-1 forms (CRF05_DF, CRF04_cpx, CRF36_cpx, CRF14_BG). These studies also revealed an unequal subtype B distribution among the different risk groups, probably due to the lack of introduction or low prevalence of this subtype in specific groups. Moreover, subtype B has a wide distribution among the PWID in Bulgaria and is the dominant subtype among MSM as it is expected. Subtype B was introduced relatively late among this group of the HIV-1 epidemic in the country.⁷

Interesting to note is that the main CRFs among studied HET individuals, including CRF01_AE and CRF02_AG, were found mostly in the same geographic regions where those clades were widespread in PWID, suggesting the existence of a bridge for viral transfer among these two groups.⁵

The present study provided an overview of the time-scaled phylogeny, population dynamics, and migration analysis of the main HIV-1 subtype B in HET individuals in Bulgaria. Moreover a comparison of the demographic and epidemiological characteristics between Bulgarian HET HIV-1 B infected individuals and those from the other HIV-1-positive populations was also performed.

The ML analysis performed as part of this study showed that, although sequences from other countries are sometimes intermixed with the Bulgarian isolates, some statistically significant clusters are including isolates only from Bulgaria. This phenomenon suggests that the HIV-1 subtype B epidemic in Bulgaria originated from limited numbers of introductions of viruses, followed by spread through local transmission networks, and this is also in line with the study by Salemi et al., which showed the dynamic nature of the Bulgarian epidemic characterized by viral inflow and outflow to/from other European countries.⁶

We estimated a mean evolutionary rate of HIV-1 B subtype pol gene sequences from Bulgarian HET populations of 1.7 × 10⁻³ substitution/site/year (95% HPD: 1.026 × 10⁻³–2.37 × 10⁻³). Our estimate is consistent with the order of magnitude of 10⁻³ expected for an HIV-1 pol gene. Moreover, our value is in line with previous estimates of the mean evolutionary rate for HIV-1 B subtype pol sequences provided by Bello et al., which estimated a mean value of 2.1 × 10⁻³ and by Hué et al., with a value of 2.55 × 10⁻³ (95% HPD: 1.74 × 10⁻³–3.51 × 10⁻³).^21,22

Regarding the time-scaled phylogeny, the origin of the root of the tree was traced back to the year 1980 (95% HPD: 1975–1986), which does not indicate that HIV-1 subtype B has been circulating in Bulgaria since that period (early 1980s), but that HIV-1 subtype B has been spreading globally since that time, with several subsequent introductions into Bulgaria. This molecular clock estimate is in line with that estimated by Salemi et al.,⁶ and is in line with the time frame of introduction of HIV-1 subtype B in North America.²³ Various introductions of HIV-1 B epidemic have since occurred between 1989 and 2007.

The evolutionary demography of the HIV-1 B isolates from Bulgaria showed an exponential growth of the number of infections from 1980s to 1990s followed by a constant phase up to about 2008; another exponential growth, from 2008 to approximately the year 2012, was found followed by a constant phase.

Regarding the migration analysis, a 31.6% of the observed viral gene flow occurred from Southwest to Northeast. A lower percentage (10.5%) of gene flow was also observed from Southeast to Southwest, suggesting the central role of these geographical groups as a potential drive of the HIV-1 subtype B Bulgarian epidemic. Other migration flow hypothesis has not been considered because it is not statistically significant.

Our molecular analysis reviewed a wide variety of different subtypes and recombinant forms of HIV-1 in the HET population. Moreover, the introductions of the main subtype B lineages into the HET population demonstrate the importance of this group for the maintenance and spread of the infection in Bulgaria. Furthermore, the undulating growth of the number of infections over the years indicates the presence of dynamic and uneven development of the epidemic in HET. Our study is of importance to analyze the events related to development of the epidemic and the role of the general population of HET with HIV-1 and in particular the introduction and spread of the main virus strain among these individuals.

The design of our study that included randomly selected specimens collected from HET patients diagnosed with HIV-1 may not be fully representative for all reported HET individuals in Bulgaria. The effect of these possible sampling biases in our analyses is unknown but may influence genotype introduction and distribution over time. In addition, any potential association of infection with demographic and epidemiological characteristics is based on self-reporting system of the patients and may be affected by recall or other biases.

Our previous studies highlighted the need for better control of the epidemic, and in 2008 a National program for prevention and control of HIV and sexual transmitted diseases was initiated in Bulgaria.^3,5,24 The aim of the program was to implement voluntary counseling and testing offices and clean needle provision to facilitate the coverage of the prevention among the most vulnerable groups in the country, including: PWID, MSM, imprisoned, sex workers, and their heterosexual partners. Using specific prevention interventions the public health response resulted in reducing HIV diagnoses in PWID since 2009.²⁵

In conclusion, our detailed phylogenetic analysis revealed high genetic diversity of different HIV-1 subtypes with uneven distribution across country. HIV-1 subtype B was introduced relatively early in the main group of HET and can play an essential role in maintaining the epidemic in the country. The migration pattern showed creeping spread of the infection throughout geographical regions.

Our study provides valuable molecular epidemiological information on HIV-1-infected HET in Bulgaria and highlights the importance of sustained molecular surveillance of the general populations to better understand and control the changing epidemic in Bulgaria.

Footnotes

Acknowledgment

This work was supported by a grant from the Ministry of Education and Science, Bulgaria. (contract: DN03/16.12.2016).

GenBank Accession Numbers

Bulgarian HIV-1 pol sequences used in this study have the following accession numbers: KY658482–KY658526, EF517418–EF517420, EF517422, JQ259068, EF517427, EF517433, EF517434, EF517436, KJ765562, EF517437, JQ259073, EF517446, JQ259074, JQ259076, JQ259078, EF517451–EF517453, JQ259082, EF517456, EF517459, EF517460, EF517462, EF517464, JQ259088, EF517465, EF517466, JQ259090, EF517467, EF517470, KJ765566, JQ259095, JQ259096, EF517473, JQ259101, JQ259102, KJ765567, EF517474, JQ259104, EF517481, EF517482, EF517484, EF517485, JQ259112–JQ259114, JQ259118, JQ259120, JQ259125, KJ765577, JQ259127, JQ259133, KJ765428, JQ259153–JQ259155, JQ259160, KJ765595, JQ259164, JQ259169, JQ259170, JQ259174, KJ765396, KJ765397, KJ765452, KJ765408, KJ765458, KJ765459, KJ765475, KJ765478, KJ765479, KJ765496, KJ765504, KJ765507, KJ765533, KJ765540, KJ765541, KJ765548, KJ765559.

Author Disclosure Statement

No competing financial interests exist.

References

Centers for Disease Control (CDC). Pneumocystis pneumonia—Los Angeles. MMWR Morb Mortal Wkly Rep, 1981; 30:250–252.

Beyrer

, Baral

, Walker

, et al.: The expanding epidemics of HIV-1 among men who have sex with men in low and middle income countries: Diversity and consistency. Epidemiol Rev, 2010; 32:137–151.

Alexiev

, Beshkov

, Shankar

, et al.: Detailed molecular epidemiologic characterization of HIV-1 infection in Bulgaria reveals broad diversity and evolving phylodynamics. PLoS One, 2013; 8:e59666.

Alexiev

, Shankar

, Wensing

AMJ

, et al.: Low HIV-1 transmitted drug resistance in Bulgaria against a background of high clade diversity. J Antimicrob Chemother, 2015; 70:1874–1880.

Alexiev

, Shankar

, Dimitrova

, et al.: Origin and spread of HIV-1 in persons who inject drugs in Bulgaria. Infect Genet Evol, 2016; 46:269–278.

Salemi

, Goodenow

, Montieri

, et al.: The HIV type 1 epidemic in Bulgaria involves multiple subtypes and is sustained by continuous viral inflow from West and East European countries. AIDS Res Hum Retroviruses, 2008; 24:771–779.

Stanojevic

, Alexiev

, Beshkov

, et al.: HIV1 molecular epidemiology in the Balkans—A melting pot for high genetic diversity. AIDS Rev, 2012; 14:28–36.

Angelis

, Albert

, Mamais

, et al. Global dispersal pattern of HIV-1 CRF01_AE: A genetic trace of human mobility related to heterosexual activities centralized in South-East Asia. J Infect Dis, 2015; 211:1735–1744.

Bozicevic

, Handanagic

, Lepej

, et al.: The emerging and re-emerging human immunodeficiency virus epidemics in Europe. Clin Microbiol Infect, 2013; 10:917–929.

10.

Struck

, Lawyer

, Ternes

, et al.: COMET: Adaptive context-based modeling for ultrafast HIV-1 subtype identification. Nucleic Acids Res, 2014; 42:e144.

11.

Strimmer

, von Haeseler

: Likelihood-mapping: A simple method to visualize phylogenetic content of a sequence alignment. PNAS, 1997; 94:6815–6819.

12.

Lo Presti

, Ciccozzi

, Cella

, et al.: Origin, evolution, and phylogeography of recent epidemic CHIKV strains. Infect Genet Evol, 2012; 12:392–398.

13.

Ciccozzi

, Ciccaglione

, Lo Presti

, et al.: Evolutionary dynamics of HBV-D1 genotype epidemic in Turkey. J Med Virol, 2014; 86:109–116.

14.

Price

, Dehal

, Arkin

: FastTree 2-approximately maximum likelihood trees for large alignments. PLoS One, 2010; 5:e9490.

15.

Drummond

, Suchard

, Xie

, Rambaut

: Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol, 2012; 29; 1969–1973.

16.

Drummond

, Rambaut

: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol, 2007; 7:214.

17.

Kass

, Raftery

: Bayes factors. J Am Stat Assoc, 1995; 90:773–795.

18.

Lo Presti

, Ciccozzi

, Cella

, et al.: Migration patterns of HIV-1 subtype B virus in Northern Italy. New Microbiol, 2013; 36:75–79.

19.

Campbell

: Chi-squared and Fisher-Irwin tests of two-by-two tables with small sample recommendations. Stat Med, 2007; 19:3661–3675.

20.

GraphPad. Graph Pad Prism v5.0. 2015. www.graphpad.com, accessed January 21, 2018 .

21.

Bello

, Eyer-Silva

, Couto-Fernandez

, Guimaraes

, Chequer Fernandez

, Teixeira

, Morgado

: Demographic history of HIV-1 subtypes B and F in Brazil. Infect Genet Evol, 2007; 7:263–270.

22.

Hué

, Pillay

, Clewley

, Pybus

: Genetic analysis reveals the complex structure of HIV-1 transmission within defined risk groups. Proc Natl Acad Sci U S A, 2005; 102:4425–4429.

23.

Robbins

, Lemey

, Pybus

, et al.: Characterization of early U.S. HIV-1 strains, dating of the introduction of HIV-1 into the U.S. and estimation of population history of the U.S. epidemic. J Virol, 2003; 77:6359–6366.

24.

Available at: www.mh.government.bg/media/filer_public/2015/04/17/programa-preventsiya-hiv-2008-2015.pdf, accessed January 21, 2018 .

25.

EMCDDA: Bulgaria Country Overview. 2015. Available at: www.emcdda.europa.eu/countries/bulgaria, accessed January 21, 2018 .

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.02 MB