Human immunodeficiency virus type 1 sub-subtype F1 pol sequences spanning the gene fragments that encode protease (PR) and the 5′-end of reverse transcriptase (RT) were collected from all patients for whom drug resistance testing was requested as part of routine clinical follow-up at six AIDS reference laboratories (ARL) in Antwerp, Brussels, Ghent, and Leuven (up until March 2015). These ARL perform approximately 90% of the HIV-1 drug resistance testing in Belgium. HIV-1 sub-subtype F1 pol sequences were requested from partners in other European and African countries where circulation of F1 infections has been reported (Palma et al., 2007; Paraschiv et al., 2009; Frentz et al., 2010; Lai et al., 2012; Ivanov et al., 2013; Sayan et al., 2013; Alexiev et al., 2015). Furthermore, all available HIV-1 sub-subtype F1 pol sequences were extracted from the Los Alamos HIV Sequence Database (LANL)¹ (Triques et al., 1999; Montano et al., 2005; Aulicino et al., 2007). Most viral genetic sequences were annotated with sampling time and location, which were used as calibration and discrete traits, respectively, in an analysis of spatiotemporal dispersion patterns (Supplementary Table S1). Patient socio-demographics and detailed sample information were collected from mandatory reporting and medical records at the participating ARL in Belgium (gender, year of birth, HIV-1 diagnosis date, infection stage at diagnosis, country of birth, probable country of infection, sexual risk factor for HIV-1 acquisition, viral load, CD4 and CD8 count associated to the HIV-1 pol sequence, and therapy initiation date). Infection stage at diagnosis was defined based upon laboratory data: HIV-1 confirmation results with a positive p24 test combined with a negative or undetermined blot or INNO-LIA HIV I/II (Fujirebio Europe, Ghent, Belgium) were classified as acute infections while results with a positive blot lacking a p31 band as recent infections (Verhofstede et al., 2017). In all other instances, at least when records were available, infections were classified as chronic. The study was approved by the Ethical Committee Research UZ/KU Leuven (reference S58359). This study did not require written informed consent by study participants as it only used data collected within routine clinical practice, but patients who opted out of research studies were not included. The entire dataset for this manuscript is not publicly available because it represents a dense sample of our epidemic and inappropriate use could endanger the privacy of the patients. A proportional down-sampling of the dataset (30%) is accessible via GenBank (accession numbers, MK458435-MK458523). Requests to access the entire dataset combined with well-defined project proposals can be mailed to the corresponding author who will submit it to the local ethical committee for ethical and legal clearance when not in conflict with already on-going research.

The quality of HIV-1 sequences that were generated by the Belgian ARL was assured by quality management systems that are conform to the ISO 15189 standard. In addition, the quality of all collected HIV-1 sequences was subsequently verified using the Quality Control tool hosted on the HIV sequence database LANL. Only sequences with no evidence of hyper-mutations and a maximum of 7 stop codons and/or frameshifts were retained for further analyses. Only pure sub-subtype F1 strains were retained in the study. Therefore, HIV-1 subtype classification was confirmed using the online automated subtyping tools COMET v1.0 (Struck et al., 2014) and Rega v3 (Pineda-Peña et al., 2013). In case of a non-F1 classification by COMET v1.0 or Rega v3 (7% discordant classifications), sequences were subsequently subjected to jpHMM (Schultz et al., 2009; version March 2015). If scored as F1 or F2 sub-subtype by jpHMM, recombination was ruled out using Simplot (Lole et al., 1999). Four HIV-1 subtype B reference sequences were added to the dataset, serving as an outgroup for subsequent phylogenetic analyses (GenBank accession numbers K03455, AY423387, AY173951, and AY331295). Nucleotide sequences were aligned using Muscle v3.7 (Edgar, 2004) and manually edited in MEGA 5 (Tamura et al., 2011). Nucleotide positions encoding for 43 surveillance drug resistance mutations (Bennett et al., 2009) were removed from the alignment to avoid the effect of drug-induced convergent evolution. Duplicate sequences were identified using the ElimDupes tool (LANL) and were removed when originating from the same sample. For patients with multiple sequences, only the first sequence was retained. The final alignment contained 1866 HIV-1 sub-subtype F1 sequences trimmed to the common length of 798 nucleotide positions, spanning sites 2253 to 2549 and 2661 to 3290 (as compared to the HIV-1 subtype B reference sequence HXB2, accession number K03455) (Supplementary Table S1).

The spatiotemporal origin of the Belgian HIV-1 sub-subtype F1 epidemic was investigated using a Bayesian approach. Prior to this analysis, we examined the accumulation of divergence over the sampling time span using Tempest v1.5 (Rambaut et al., 2016) to assure that sampling dates could be used as calibration. Markov chain Monte Carlo analyses were performed using BEAST v1.8.4 in conjunction with BEAGLE to improve the computational performance (Ayres et al., 2012; Suchard et al., 2018). We specified a GTR nucleotide substitution model with discrete gamma distributed rate variation among sites. Among-lineage variation in evolutionary rates was modeled using a lognormal relaxed molecular clock and a flexible non-parametric Bayesian Skygrid was specified as a coalescent prior. The analyses were run in triplicate for 100 million states and trees were sampled every 10,000th state. Convergence and mixing properties of the Markov chain Monte Carlo chains were assessed with Tracer v1.6 (Rambaut et al., 2018). Results of several chains were combined after removal of the burn-in, and a MCC tree was summarized using TreeAnnotator within the BEAST package. The timed migration process was reconstructed using a model that allows for non-reversible migration rates. Trees were visualized with FigTree v1.4.3 and iTOL v4.2.3 (Letunic and Bork, 2016).

Phylodynamic analyses were performed on an infection cluster that mainly contained sequences from our Belgian cohort. A birth-death model was applied to the Belgian sequences included in this infection cluster (191 sequences) in order to quantify transmission dynamics through time. The analyses were performed in BEAST v2 (Bouckaert et al., 2014) using the BDSKY add-on (Stadler et al., 2013). The molecular evolution and clock model were chosen as in the phylogeographic analysis. For the birth-death-skyline parameters, we used a LogNorm(0,1) prior for the effective reproductive number, a LogNorm(0,1) prior for the becoming-non-infectious rate in units per year (i.e., the inverse of the time duration of being infectious in unit of years), and a Unif(0,infinity) prior for the origin of the epidemic. The dimension of the effective reproductive number and the becoming-non-infectious rate was set to 10, meaning the parameters were estimated for 10 equally spaced intervals between the origin and the most recent sample. The sampling proportion was fixed to 0.65 (and 0.8 in a second analysis) for the time between the first and last sample, and 0 prior to the first sample.

Socio-demographic, virological and clinical variables were reported as absolute numbers, percentages and medians with their respective interquartile ranges (IQR), where appropriate. Their association with a Belgian F1 outbreak was evaluated using univariate and multivariate logistic regression in the statistical R software v3.5.1, with statistical significance set at p-value < 0.05. For these analyses, only patients with available information for a specific variable were included.

Up to February 2015, 297 HIV-1 sub-subtype F1 infections were newly diagnosed in the participating Belgian ARL (Figure 1), of which 137 were analyzed in Antwerp (46%), 95 in Brussels (32%), 38 in Leuven (13%), and 27 in Ghent (9%). The median age of people diagnosed with an HIV-1 sub-subtype F1 infection was 34 years (IQR 29–41) and 77% were men. Over the course of 15 years, the number of HIV-1 sub-subtype F1 diagnoses increased from 3 in March 2001 – February 2003 to 88 in March 2013 – February 2015 (Figure 1A). Sub-subtype F1 infections were initially mainly diagnosed in women (65%), but after February 2005 men became the most affected gender (87%). Of the sub-subtype F1 infections for which the respective information was registered, the majority were diagnosed in people born in Belgium (59%, 141/238) and in people who declared that they probably acquired their infection in Belgium (82%, 138/168). The sub-subtype F1 infections were found in homosexual (51%, 117/231), heterosexual (35%, 80/231), and bisexual risk (15%, 34/231) groups (Figure 1B).

Bayesian phylogenetic analyses were used to estimate the geographical and temporal origin of the HIV-1 sub-subtype F1 strains in our national cohort and to investigate viral dissemination from Belgium to other countries. The final dataset included 1866 sequences that were mainly sampled in Europe (75%) and Brazil (20%) (Supplementary Table S1). We found sufficient temporal signal to pursue phylodynamic analysis and used the resulting empirical distribution of the dated phylogenetic trees in a subsequently phylogeographic analysis.

The phylogeographic analysis revealed the divergence of sub-subtype F1 into two well-supported clades for which (i) one clade included mainly Brazilian and Portuguese strains together with a Belgian subclade, while (ii) the other predominantly Belgian, Romanian and Portuguese sequences (upper and lower clade in Figure 2, Supplementary Figure S1). Their divergence dates back 50 years ago (1968, 95% higher posterior density (HPD) interval of 1963–1973), with an estimated start of local dissemination of HIV-1 F1 in Brazil around 1973 (95% HPD 1969–1975) and subsequently in Portugal around 1979 (95% HPD 1976–1981) (Figure 2 upper clade). In this upper clade (i), one large Belgian cluster, subsequently referred to as the sub-subtype F1 Belgian outbreak, is embedded within a Western-European clade that in turn is nested within the larger Brazilian clade. The Belgian clade contains four sequences sampled outside Belgium (one from Brazil, one from Luxembourg, and two from Germany), while the Western-European clade consists of the Belgian clade together with three basal sequences from Switzerland, one large well-supported Spanish clade and one Romanian sequence that branches outside the Spanish and Belgian clade. The founder of the HIV-1 sub-subtype F1 strains that are currently circulating in Western-Europe, probably migrated from Brazil to Europe between 1990 and 1999 (95% HPD 1987–2002), whereas the origin of the Belgian and Spanish outbreaks were dated back to 2002 (95% HPD 2000–2004) and 2006 (95% HPD 2004–2008), respectively. Brazil is estimated as the most likely geographical origin of the Western-European clade because the Western-European clade was embedded within the large Brazilian sub-subtype F1 diversity and because one strain that was sampled in Brazil and three strains that were sampled from Brazilians residing in Belgium were also positioned within the Belgian clade. Although Belgium was assigned as ancestral location state of this Western-European clade, the long interconnecting branch suggests that current sampling is insufficient to provide conclusive results as to the specific migration event(s) that finally led to the two large outbreaks in Western Europe (Belgium and Spain).

The remaining sub-subtype F1 sequences from our Belgian cohort represent single lineages or small clusters that are dispersed across the tree and are intermixed among viral lineages from Africa, Latin-America and Europe (Figure 2, 3). Several other strains from Sub-Saharan migrants residing in Belgium but who originated from the Democratic Republic of Congo (DRC), Congo, and Angola occupied the most basal positions of the tree (Figure 2 lower clade, Figure 3). They cluster with African lineages that probably migrated to Europe between 1972 and 1977 (95% HPD 1969–1980) and that were at the origin of several larger sub-subtype F1 outbreaks in Europe, more specifically in Portugal (1979, 95% HPD 1976–1982), Romania (1981, 95% HPD 1978–1983; 2001, 95% HPD 1996–2006; and 2002, 95% HPD 1998–2006), and Bulgaria (1982, HPD 1979–1985). As, we defined a location trait according to country of sampling in our analysis, the over-representation of Belgian sequences – probably to a large extent reflecting African diversity – led Belgium to be inappropriately assigned as the ancestor location of this lower clade that contained many African, Romanian, Portuguese, and other European lineages (Figure 2 lower clade).

As one or only a few sub-subtype F1 lineages have successfully established local transmissions among 191 of 297 people diagnosed with a sub-subtype F1 infection in Belgium (64%), we used univariate and multivariate logistic regression to identify socio-epidemiological variables that were associated with this successful sub-subtype F1 outbreak in Belgium (Table 1). HIV-1 patients who were diagnosed with a sub-subtype F1 infection at the AIDS reference laboratories in Antwerp (OR 6.80, 95% CI 3.86–11.99), Leuven (OR 4.22, 95% CI 1.59–11.16) or more broadly in the region of Flanders (OR 7.28, 95% CI 4.23–12.53) had significantly higher odds of belonging to this F1 outbreak, compared to people who were diagnosed with a sub-subtype F1 infection elsewhere in Belgium. In addition, many other variables of the sub-subtype F1 infections diagnosed in Belgium displayed significantly higher odds of being part of the outbreak, of which the most informative were male gender (OR 52.42, 95% CI 19.90–138.08), infection in MSM (combining homosexuals and bisexuals) (OR 34.8, 95% CI 16.0–75.6), probable infection in Belgium (OR 26.67, 95% CI 9.59–74.15), born in Belgium (OR 12.10, 95% CI 6.35-23.08), recent infection stage at diagnosis (OR 7.19, 95% CI 2.88–17.95) and higher viral load (OR 2.36, 95% CI 1.78–3.13). In multivariate analysis, sub-subtype F1 infections in MSM and in people born in Belgium remained to display a higher likelihood for belonging to the sub-subtype F1 Belgian outbreak [OR 7.82, 95% CI 3.80–16.11 and OR 2.01 (1.06–3.78), respectively]. The variables sampling region, ARL in Belgium, gender, country of birth and sexual risk group were visualized for all Belgian sequences in a circular representation of the phylogenetic tree (Figure 3).

Since the onset of the sub-subtype F1 outbreak in Belgium (2002, 95% HPD 2000–2004), we covered a time period where antiretroviral therapy was initially still deferred in asymptomatic patients with CD4 count above 350 cells/mm³ (European AIDS Clinical Society, 2005). This general practice gradually switched to the consideration of antiretroviral therapy for asymptomatic patients with CD4 counts between 350 and 500 cells/mm³ (European AIDS Clinical Society, 2007) and above 500 cells/mm³ (European AIDS Clinical Society, 2013). Given the benefits of antiretroviral therapy in reducing infectivity, we wanted to quantify transmission dynamics through time using a birth-death model and monitor the trends in therapy uptake in our specific setting. Therapy and diagnosis information was available for 147 HIV-1 patients whose viral strains belonged to the sub-subtype F1 outbreak in Belgium (77%, 147/191). The delay in therapy initiation decreased from a median of 1011 days (IQR 368–1683) or almost 3 years for patients diagnosed in 2005 and 2006 to a median of 61 days (IQR 26–208) or about 2 months for patients diagnosed in 2013 and 2014 (Figure 4). In the birth-death model, assuming a constant 65% sampling proportion, the number of expected secondary infections per infected person (effective reproductive number R_e) was estimated to be a median of 3 in 2005 and gradually declined to a median of 1.5 in 2011, after which a substantial drop to below 1 was observed (Figure 5). During this time window, the median infectious time also declined and the number of expected secondary infections per infected person and per year (transmission rate) decreased from a median of 1 to well below 1. Sensitivity analysis using a constant sampling proportion of 80%, as the diagnosis and sequencing rate is expected to be higher in MSM, showed similar results and suggests that the estimates are relatively robust to prior specification on the sampling proportion.

Our and other phylogenetic analyses confirmed that the most basal F1 sequences were isolated from people originating from Central African countries and that the major migration events out of Africa occurred in the 1970s and early 1980s seeding some larger F1 epidemics over the world (Mehta et al., 2011; Bello et al., 2012; Lai et al., 2014). During the Angolan civil war that started in 1975 and lasted for more than 16 years, an extensive migration of troops, settlers and refugees was observed between Angola, DRC, Portugal, Brazil, Cuba, Russia, and several Eastern-European countries and this is considered as an important factor in the global dissemination of F1. In Brazil, sub-subtype F1 infections have reached prevalence rates of more than 7% and the strains are currently circulating among various risk groups (Brindeiro et al., 2003). It is the most probable epicenter of F1-related epidemics within the Southern Cone and one of the sources of the Portuguese sub-subtype F1 epidemic (Aulicino et al., 2007; Bello et al., 2007, 2012; Carvalho et al., 2015). Unrelated to the Brazilian epidemic, sub-subtype F1 was initially introduced and spread among adult heterosexuals in Romania. In this Eastern-European country, it subsequently reached a very high prevalence rate of more than 70% due to several major outbreaks among children in orphanages and people using intravenous drugs (Dumitrescu et al., 1994; Bandea et al., 1995; Apetrei et al., 1998; Op De Coul et al., 2000; Guimarães et al., 2009; Paraschiv et al., 2009; Mehta et al., 2011; Mbisa et al., 2012; Niculescu et al., 2014). Both Brazil and Romania were at the origin of the F1 epidemic in Italy, where Brazilian strains contributed to many more introduction events characterized with a more extensive spread (Lai et al., 2014). Since then, the prevalence of sub-subtype F1 increased gradually among Italian descendants and mainly among self-declared heterosexual men, although phylogenetic analyses suggested that some of these self-declared heterosexual men might have been infected through homosexual contact (Bracciale et al., 2009;Lai et al., 2012, 2014).

Although our and a previous study indicated that another sub-subtype F1 lineage migrated from Brazil to Europe during the 1990s, this lineage obtained only more recently epidemic proportions in two countries, respectively, about 16 years ago in Belgium and 12 years ago in Spain (Paraskevis et al., 2015). In both instances, F1 was introduced in local MSM networks, a factor probably important in its subsequent spread and circulation at high rates (Thomson et al., 2012). Our analyses did not indicate any intermixing between both epidemics and only sporadic and putative exports from Belgium to Luxembourg (one) and Germany (two) were detected. This is in agreement with another study that traced the HIV-1 subtype B mobility in Europe and that characterized Belgium primarily as a sink of HIV-1 lineages rather than a source (Paraskevis et al., 2009). Although, we confirmed previous findings on the global spread of HIV-1 sub-subtype F1, our sampling approach might not have been sensitive enough to capture the specific migration events that resulted into these two large outbreaks or into extensive onward transmission to other countries. For our analyses, we collected HIV-1 pol sequences from countries with a known high burden of F1 infections or a potential role in the global spread of sub-subtype F1. In this dataset, Belgium, Romania, Portugal, Luxembourg, and Bulgaria were the best-represented countries as they contributed a dense sampling, whereas Italy provided a down-sampling proportional to their local sub-subtype F1 diversity. However, other countries might have been undersampled as they were represented by sequences that were received from single cohorts or by sequences that were collected from public databases. Such sample biases call for caution in interpreting the results of phylogeographic estimates.

The first HIV-1 sub-subtype F1 infections within the local outbreak cluster were diagnosed in 2005. It concerned one chronic infection in Leuven and one infection with unknown status in Antwerp, suggesting a time of about 3 years between the estimated introduction in Belgium (2002, 95% HPD 2000–2004) and the first case reported by the Belgian ARL. Although there is universal access to HIV-1 health care in Belgium and we included all patients for whom genotypic drug resistance testing indicated a sub-subtype F1 infection, whether at baseline or at therapeutic failure, our sequencing efforts in routine practice only cover approximately 60% of all newly HIV-1 diagnoses registered in the national database. However, we believe that we do capture the majority of patients who are contributing to the local spread, as sub-subtype F1 is mainly circulating among MSM, a group who is well engaged in medical care (Van Beckhoven et al., 2015). Any HIV-1 patient who seeks medical care in Belgium will receive a confirmation of his/her infection and will therefore be registered in the national database. Therefore, we believe that a large proportion of the people living with HIV-1 in Belgium who do not receive a baseline resistance test at diagnosis are people with an infection that was already acquired and diagnosed abroad and who are already on successful antiretroviral treatment upon arrival in Belgium, or these people are only temporarily in the country and do not enter the HIV-1 health care system in Belgium.

With the dense sampling of our local epidemic at hand, we aimed to estimate the effective reproductive number R_e, the transmission (birth) rates and the time of being infectious in our viral population using birth-death modeling in order to gain insights into the dynamics of the outbreak. The initial high R_e-value, high transmission rate and long times of being infectious were in agreement with the explosive increase of sub-subtype F1 infections in our national cohort. However, the mean R_e-value declined substantially from 2012 onward and this event was associated with lower transmission rates and shorter durations of being infectious. These observations may be explained by the policy changes around that time when therapy eligibility thresholds for asymptomatic HIV-1 patients became less strict in Belgium. Indeed, medical records showed a remarkable dent in the median treatment delay among people belonging to our outbreak around 2009–2010. Although this relationship does not demonstrate causality, it is supported by studies where viral suppression by antiretroviral therapy has shown to be effective in preventing HIV-1 transmission among serodiscordant heterosexual and MSM couples (Cohen et al., 2011; Rodger et al., 2016; Bavinton et al., 2018). However, the relevant changes in transmission dynamics were not yet reflected in the number of new sub-subtype F1 diagnoses up to the start of 2015. This could suggest that the effects were too recent to be fully captured by our surveillance data of new diagnoses or that the conditions were not yet fully optimal to achieve the full potential of treatment as prevention in this outbreak (Okano et al., 2016). Although a substantial proportion of our outbreak was diagnosed at an acute and recent stage of infection, we should aim to increase testing and diagnosis rates among people at risk of acquiring HIV-1 and treat them without any delay. In this respect, the INSIGHT START study resolved the last concerns for immediate initiation of antiretroviral therapy in asymptomatic HIV-1 infections as this study showed that this strategy was superior to deferral until the CD4 count declined to 350 cells/mm³ and that it was not associated with an increased rate of adverse events (INSIGHT START Study Group et al., 2015).

In conclusion, this study revealed that sub-subtype F1 strains were introduced in Belgium around 2002 and that during a decade they spread exponentially among MSM who were mainly from Belgian descent, residing in Flanders and diagnosed at a recent stage of infection. Our birth-death phylodynamic estimates are in line with the expectations of the effectiveness of treatment as prevention because we identify a decline in the onward transmission of sub-subtype F1 strains that coincides with the implementation of earlier initiation of antiretroviral therapy in our cohort. Further follow-up is planned, especially within the context of the roll-out of antiretroviral therapy in early asymptomatic HIV-1 infections and of the implementation of pre-exposure prophylaxis which is reimbursed to people who are at increased risk for HIV-1 infection since June 2017.