TRACING THE ORIGIN AND DYNAMICS OF THE HIV-1 EPIDEMIC IN SERBIA

Since the first report of HIV infection in Serbia in 1985, the HIV-1 epidemic was very dynamic, changing the pattern in subtype distribution and prevailing transmission routes. To better understand the origin and epidemiological dynamics of HIV-1, we analyzed 266 (pol) sequences from Serbian patients diagnosed over a period of 14 years. Subtype distribution in Serbia is still marked by a prevailing subtype B genetic form. The transmission pattern, however, has changed from being intravenous drug user (IVDU) driven to predominantly sexual transmission. The estimated time of initial founder strain introduction of sequences from Serbian IVDUs and MSM (men who have sex with men) is similar and dates back to the early 1980s, while introduction of subtype C occurred much more recently. Key words, Serbian epidemic, subtype distribution, HIV-1 transmission, Bayesian coalescent analyses


INTRODUCTION
The high genetic variability and rapid evolution of HIV-1 are important factors in its worldwide spread.HIV genetic heterogeneity originates from the high mutation and recombination rates of the reverse transcriptase enzyme combined with a high turnover rate.This results in genetically diverse populations of viral species in each infected individual, termed 'quasispecies' .Most HIV-1 infections globally are caused by group M viruses, with groups O, N and P causing a small minority of infections in central Africa (Peeters and Sharp, 2000).Currently, within the group M, there are nine recognized subtypes of HIV-1, a great variety of circulating recombinant forms and unique recombinant forms (McCutchan, 2000;Robertson et al., 2000;Hemelaar et al., 2012).While subtype B predominates in North America, Australia, Western and Central Europe (Peeters et al., 2000), it is estimated that more than 55% of worldwide HIV-1 infections are caused by HIV-1 subtype C.
The HIV-1 epidemic in Serbia is mainly driven by the subtype B genetic form, even though other and novel subtypes have been identified.Among non-B subtypes, subtype G was found to be the most prevalent one in previous studies, accounting for almost one third of all non-B infections.(Stanojevic et al., 2012;Siljic et al., in press).Among European countries, Portugal is the only one where this subtype is found in an important proportion of HIVinfected patients (Esteves et al., 2003;Abecasis et al., 2013).Other non-B subtypes and circulating recombinant forms (C, A, F, CRFs) were found to be less represented (Stanojevic et al., 2002;Stanojevic et al., 2012).
Serbia has around 3 000 registered cases of individuals living with HIV-1 (PLWH), being a country with a low prevalence of HIV infection (http://www.batut.org.rs/download/aktuelno/ Epidemija %20 HIV% 20infekcije%20u%20Srbiji%202012.pdf).The first case of HIV infection officially reported in Serbia dates back to 1985.So far, epidemiological and existing molecular data indicate that the local HIV-1 epidemic began with subtype B, among intravenous drug users (IVDU), and this transmission route was the most prevalent one over subsequent years (Stanojevic et al., 2002).However, Serbia has recently experienced a critical change in its HIV-1/AIDS epidemic, with significant and rapid increase in the rates of HIV-1 transmission through sexual contact (http:// www.batut.org.rs/download/aktuelno/ Epidemija %20HIV% 20infekcije%20u%20Srbiji%202012.pdf).Particularly concerning is the increasing trend of infection among men who have sex with men (MSM).
Reconstructing the history of the introduction of the founder strain is crucial in understanding the actual characteristics of the Serbian HIV-1 epidemic.According to our previous findings, the epidemic spread of HIV in Serbia within the local network and outside the firstly implicated IVDU community probably started at the beginning of the 1990s, and the estimated origin of subtype G epidemics in Serbia has also been dated to the beginning of the 1990s (Siljic et al., in press).However, estimates about the initial introduction of HIV-1 B in Serbia as well as other genotype viruses that are commonly found in Serbia still remain uncertain.To fill this gap, we reconstructed the time to the most recent common ancestor (tMRCA) of the predominant HIV-1B strains detected in PLWH with different risk for HIV infection, as well as subtype C isolates from Serbia, as the second most prevalent non-B subtype in the country (Siljic et al., in press).

Study population
Blood samples used in the study were collected from 1997 to 2011 from 266 HIV-1 seropositive patients attending the HIV/AIDS Center, Institute of Tropical and Infectious Diseases, Belgrade, Serbia.Informed consent was obtained from each patient.Consenting individuals were above 18 years of age, either antiretroviral naïve or treated, with viremia above 1 000 copies per ml and presenting with different risks for acquiring HIV infection.Transmission risk was categorized as men who have sex with men (MSM), heterosexual, intravenous drug use (IVDU), transfusion, vertical transmission or unknown.For patients reporting IVDU in addition to another risk, the former was considered as main risk for acquiring HIV infection.The patient epidemiological, demographical and clinical information included age, gender, place of birth, route of infection, viral load and CD4 cell count.

RNA extraction, amplification and DNA sequencing
Plasma was separated within six hours after collection and frozen at -80 0 C for further analysis.Viral RNA was isolated from 280 μL of plasma using a QIAamp viral RNA kit (Qiagen Inc., Germany) according to the manufacturer's instructions.The 1.6-kb pol region (entire protease and 256 codons of the reverse transcriptase coding gene) were reverse-transcribed and amplified with the One-Step RNA PCR Kit (Qiagen, Hilden, Germany) using previously described primers and methods.In the second round of PCR, amplification was performed with a Taq PCR Core Kit (Qiagen, Hilden, Germany).Nested PCR products were purified using a Qiagen PCR purification kit (Qiagen Inc.) according to the manufacturer's protocol.Purified DNA was sequenced directly in both directions using six sequencing primers and the ABI BigDye Terminator v.3.1 cycle sequencing ready reaction kit (Applied Biosystems, Foster City, California, USA) and processed with an automated ABI 310 Genetic Analyzer (Applied Biosystems) (Snoeck et al., 2005).

Sequence assembly, quality control and phylogenetic analysis
The sequences were assembled using SeqScape HIV-1 Genotyping System Software v 2.5 (Applied Biosys-tem, Foster City, CA, USA).All assembled sequences were aligned in ClustalW as implemented in MEGA 5.1 software (http://www.megasoftware.net/), then visually inspected and manually edited when necessary.All positions that contained alignment gaps were removed.jModelTest program was used to select the best-fit model for nucleotide substitution, resulting in the GTR+I+Γ model in the dataset (Posada et al., 2008).A maximum-likelihood tree was constructed under the selected model, as implemented in PAUP v4.0beta10 (Swofford, 1998).The topology of trees was tested by bootstrap analysis with 1 000 replicates.The tree was edited and displayed using FigTree v1.3.1.

HIV-1 subtyping
Subtype assignment was performed on all 266 sequences obtained in the study using the REGA HIV-1 subtyping tool (http://jose.med.kuleuven.be/genotypetool /html/subtypinghiv.html).Furthermore, 250/266 Serbian sequences of sufficient length (fulllength protease and minimally 250 codons of reverse transcriptase) were used in ML phylogenetic analysis together with the reference sequences downloaded from the HIV Los Alamos database (www.hiv.lanl.gov/content/index) for additional subtype assignment and general phylogenetic insight into the dataset.

PhyloPart phylogeny of Serbian MSM subtype B HIV-1sequence
A PhyloPart transmission cluster analysis was performed on a maximum-likelihood phylogeny tree with 121pol MSM sequences of sufficient length, reconstructed for this purpose under the selected best-fit model GTR+I+Γ as implemented in PAUP.Presence of transmission clusters was assessed under the following parameters: threshold 0.05 and sample distance limit 6 (Prosperi et al., 2011).

Evolutionary analysis
Estimates of the time of the most recent common ancestor (tMRCA) for the sequences from IVDUs,  (Shapiro et al., 2006;Abecasis et al., 2009).MCMC chains were run for 10 million generations and sampled every 1 000 steps.Bayesian MCMC output was analyzed using TRACER v1.5 (http://tree.bio.ed.ac.uk/ software/ tracer/), with uncertainty in parameter estimates reflected by the 95% highest probability density (HPD) values.The Effective Sample Size (ESS) values for estimates were more than 100.The trees were consolidated into a target tree using the TreeAnnotator program and scanned using the FigTree program.

Statistical analyses
The chi-square test and Fisher's exact test were used for statistical analyses.P values <0.05 were considered statistically significant and P values <0.01 were considered highly statistically significant.

RESULTS
The general epidemiologic data of the study population are shown in Table 1.The present study included 78% (208/266) of men, of whom the majority, 59% (123/208), reported MSM as the route of HIV-1 transmission.In the overall populations, almost half of the patients, 46.2% (123/266), reported an MSM risk for acquiring HIV infection, followed by heterosexual contact in 31.9%(80/266) and intravenous drug use in 6.8% (18/266).Blood transfusion and vertical transmission were the reported risk for infection in a very small proportion of all patients, 4.5% (12/266) and 1.1% (3/266), respectively.
Both the Rega subtyping tool and phylogenetic analysis gave congruent results in terms of subtype assignment (Fig. 1).The major circulating subtype in Serbia is subtype B found in 91.4% (243/266) of the collected samples.Among non-B subtypes, subtypes G and C were the most prevalent, found in 2.3% (6/266), each.Subtype A and subtype F were found in a very small proportion, 1.3% (3/266) and 0.4 % (1/266), respectively.Circulating recombinant forms, CRF01_AE and CRF02_AG, were found in 3.2% (7/266).
To understand better the temporal trends of the epidemiology of HIV-1 infection in Serbia, we analyzed trends in the prevalence of HIV-1 subtypes and trends in the major risk groups across two equal time periods: 1997-2004 (136 sequences sampled) in comparison to 2005-2011 (130 sequences sampled).During the first half of the present survey (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004), intravenous drug use was a significantly more frequently reported transmission risk (p<0.001)than in the second half (2005)(2006)(2007)(2008)(2009)(2010)(2011).In contrast, MSM was the significantly more prevalent reported transmission risk (p<0.0001) in the second half of the survey.Regarding the changing pattern in HIV-1 subtype  distribution, we found a highly significant increase in the prevalence of subtype B genetic form in the period from 2005 till the end of the present research (p<0.01),while in the same period neither subtype G nor subtype C were identified in any of the samples collected in the present study.
To determine the epidemic origins of HIV strains among Serbian IVDUs, we estimated the tMRCA by using the relaxed molecular clock model based on the nucleotide sequences of 1.6-kbpol gene.The tMRCA inferred for local IVDU sequences was in 1983 (95% Higher Posterior Density HPD: 1974-1993) (Fig. 2).Evolutionary analyses of pol genomic regions indicate that the age of the most recent common ancestor for the HIV-1 subtype B sequences circulating among Serbian MSM dates back to 1984 (95% HPD: 1979-1989) (Fig. 3).The estimated temporal origin for the local subtype C lineage was much more recent, 1993recent, (95% HPD: 1987recent, -1999)).

DISCUSSION AND RESULTS
Since the first reported HIV-1-infected case in Serbia, subtype B has remained the predominant one.However, with changes in the distribution of non-B subtypes over the years, new non-B subtypes and increased genetic diversity among them has emerged.This situation is similar to western and central European countries with the highest prevalence of subtype B observed in Poland (96.2%) and Slovenia (93.6%), while the lowest proportion of subtype B was found in Israel (27.9%) and Portugal (39.2%) (Paraskevis et al., 2009;Abecasis et al., 2013).The most common non-B subtypes in Serbia were found to be subtypes G and C. Worldwide, HIV-1 subtype C is currently the most prevalent subtype, accounting for approximately 55% of all infections.Among the countries of the WHO European region, subtype C is the most prevalent in Israel, where it accounts for 58% of all HIV infections, while subtype G is the most prevalent in Portugal (34.8%),where it is well established in the native population (Abecasis et al., 2013).Previous study concerning the molecular epidemiology of HIV-1 in Serbia found that non-B subtypes are more frequently seen in the heterosexual population, while subtype B is more common in MSM (Siljic et al., in press).
Here we describe the changing epidemiology of HIV-1 disease among main transmission groups.In Serbia, as in other countries, the HIV epidemic is evolving with changing patterns in the HIV subtypes present and prevailing transmission risk.In recent years, the majority of men with newly diagnosed HIV infection in Serbia reported MSM contact as the most probable risk for HIV infection.According to epidemiological data, MSM contact currently represents the driving force of the HIV epidemic in Serbia, followed by heterosexual contact, while infections among injecting drug users are at low levels.
The availability of greater numbers of sampled viral sequences combined with the increasing power of phylogenic techniques makes it possible to retrieve valuable and unique information about the course of viral epidemics from molecular data (Lemey et al., 2003;Robbins et al., 2003;Travers et al., 2004;Deng et al., 2008).Since primo-infection, the genetic distances from the founder strain of the HIV-1 quasispecies within an individual increase in a linear fashion over time (Shankarappa et al., 1999;Bello et al., 2006).At the population level, the mean rate of divergence of HIV-1 has been used to estimate the date of origin of the HIV-1 and HIV-2 pandemics (Lemey et al., 2003;Robbins et al., 2003), as well as the origin of several local epidemics (Abebe et al., 2001;Lukashov et al., 2002;Paraschiv et al., 2012;Yebra et al., 2013).In the present study, we used phylogenetic analyses and a Bayesian coalescent-based framework to investigate the origin and to estimate with more precision the initial introduction of the HIV-1 subtype B in Serbia.We based this analysis on datasets of pol gene sequences sampled over a period of 14 years.Molecular clock analysis has dated the initial introduction of both the founder strain of IVDUs and MSM to the early 1980s.Our findings suggest that MSM and IVDU populations might have played a major role in the introduction and initial dissemination of HIV-1 in Serbia.The initial introduction during the early phase of the HIV-1 epidemic in Serbia was followed by subsequent regional spread among the MSM population that gave rise to local transmission chains.The estimated origin of the subtype C epidemics in Serbia is much more recent, however, being similar to the time found for subtype G in a previous study (Siljic et al., in press).
Studies of the evolutionary history of HIV enable us to explore the past and to estimate the age of an epidemic.Phylogenetic analysis has successfully been used to investigate HIV-1 transmission and epidemiological linkage; however, the direction and timing of HIV-1 transmission are far more difficult to assess (Rachinger et al., 2011).Standard phylogenetic reconstruction only establishes evolutionary relatedness but not the evolutionary direction.The integration of molecular clock models in phylogenetic inference has enabled the reconstruction of rooted, time-measured phylogenies (Drummond et al., 2002;Rachinger et al., 2011).Estimating the date of the most recent common ancestor of viral strains, including dating the origin of HIV-1 and HIV-2, has proven to be possible with molecular data for wellcharacterized sequences with known sampling dates (Korber et al., 2000;Salemi et al., 2001;Lemey et al., 2003;Travers et al., 2004).
Worldwide, several studies estimating the date of HIV-1 origin have been done so far (Korber et al., 2000;Robbins et al., 2003;Hué et al., 2005;Gilbert et al., 2007).The estimated coalescence time of the simian ancestor of HIV-1 group M and its closest related cpz strains occurred around the end of the XVII century, while the calculated origin of HIV-1 group M radiation was in the 1930s (Korber et al., 2000;Salemi et al., 2001).Subsequently, this founder strain diversified into the various subtypes that we see worldwide today.The first subtype B migration event out of Africa is estimated to have occurred around 1966 in Haiti (Gilbert et al., 2007).The virus then moved at the end of the 1960s to the USA (Lukashov and Goudsmit, 2002).
Our data support increased public health interventions targeting MSM.Founder strains of both IVDUs and MSM were imported into Serbia in the similar period, 1983 and 1984, respectively.However, there has been a highly significant increase in the incidence of HIV infection in MSM between 2004 and 2011, reflected in the increasing number of local transmission chains, especially among middle-aged men living in Belgrade.A significant increase in the subtype B genetic form as well as absence of further spread of subtype G and subtype C epidemics in the second period of our research are in positive association with the changing epidemiology of the HIV-1 epidemic among different risk groups.

Fig. 2 .
Fig. 2. Estimated temporal origin of Serbian IVDUs sequences.Reference sequences sampled across Europe America with knowing collection date are used as an outgroup.The scale is in years.
Serbian MSM subtype B HIV-1sequences revealed the existence of three major transmission clusters.Transmission clusters included 38, 8 and 5 viral sequences.The median age of the patients in the clusters was 33.The majority of viral sequences in all three major transmission clusters is from patients living in Belgrade, 68.2% (35/51).The most recent common ancestors (tMR-CA) of analyzed Serbian MSM transmission clusters were in 1989 (95% HPD: 1980-1998) for the most expanded transmission cluster; 1997 (95% HPD: 1990-2004) for transmission cluster consisting of 8 viral sequences, and 2001 (95% HPD: 1993-2009) for transmission cluster consisting of 5 viral sequences.