Front. Microbiol. Frontiers in Microbiology Front. Microbiol. 1664-302X Frontiers Media S.A. 10.3389/fmicb.2022.821006 Microbiology Original Research Molecular Epidemiology and Transmission Dynamics of the HIV-1 Epidemic in Ethiopia: Epidemic Decline Coincided With Behavioral Interventions Before ART Scale-Up Arimide Dawit Assefa 1 2 Esquivel-Gómez Luis Roger 3 Kebede Yenew 4 Sasinovich Sviataslau 1 Balcha Taye 1 Björkman Per 1 Kühnert Denise 3 Medstrand Patrik 1 * 1Department of Translational Medicine, Lund University, Malmo, Sweden 2TB/HIV Department, Ethiopian Public Health Institute, Addis Ababa, Ethiopia 3Transmission, Infection, Diversification and Evolution Group, Max-Planck Institute for the Science of Human History, Jena, Germany 4Africa Centre for Disease Prevention and Control, Africa Union Commission, Addis Ababa, Ethiopia

Edited by: Kok Keng Tee, University of Malaya, Malaysia

Reviewed by: Teiichiro Shiino, National Center for Global Health and Medicine, Japan; Edson Delatorre, Federal University of Espirito Santo, Brazil; Kledoaldo Lima; Hospital das Clínicas, Universidade Federal de Pernambuco, Brazil

 *Correspondence: Patrik Medstrand, patrik.medstrand@med.lu.se

This article was submitted to Virology, a section of the journal Frontiers in Microbiology

 25 02 2022 2022 13 821006 23 11 2021 24 01 2022 Copyright © 2022 Arimide, Esquivel-Gómez, Kebede, Sasinovich, Balcha, Björkman, Kühnert and Medstrand. 2022 Arimide, Esquivel-Gómez, Kebede, Sasinovich, Balcha, Björkman, Kühnert and Medstrand

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

 Background

Ethiopia is one of the sub-Saharan countries hit hard by the HIV epidemic. Previous studies have shown that subtype C dominates the Ethiopian HIV-1 epidemic, but the evolutionary and temporal dynamics of HIV-1 in Ethiopia have not been closely scrutinized. Understanding the evolutionary and epidemiological pattern of HIV is vital to monitor the spread, evaluate and implement HIV prevention strategies.

 Methods

We analyzed 1,276 Ethiopian HIV-1 subtype C polymerase (pol sequences), including 144 newly generated sequences, collected from different parts of the country from 1986 to 2017. We employed state-of-art maximum likelihood and Bayesian phylodynamic analyses to comprehensively describe the evolutionary dynamics of the HIV-1 epidemic in Ethiopia. We used Bayesian phylodynamic models to estimate the dynamics of the effective population size (Ne) and reproductive numbers (Re) through time for the HIV epidemic in Ethiopia.

 Results

Our analysis revealed that the Ethiopian HIV-1 epidemic originated from two independent introductions at the beginning of the 1970s and 1980s from eastern and southern African countries, respectively, followed by epidemic growth reaching its maximum in the early 1990s. We identified three large clusters with a majority of Ethiopian sequences. Phylodynamic analyses revealed that all three clusters were characterized by high transmission rates during the early epidemic, followed by a decline in HIV-1 transmissions after 1990. Re was high (4–6) during the earlier time of the epidemic but dropped significantly and remained low (Re < 1) after the mid-1990. Similarly, with an expected shift in time, the effective population size (Ne) steadily increased until the beginning of 2000, followed by a decline and stabilization until recent years. The phylodynamic analyses corroborated the modeled UNAIDS incidence and prevalence estimates.

 Conclusion

The rapid decline in the HIV epidemic took place a decade before introducing antiretroviral therapy in Ethiopia and coincided with early behavioral, preventive, and awareness interventions implemented in the country. Our findings highlight the importance of behavioral interventions and antiretroviral therapy scale-up to halt and maintain HIV transmissions at low levels (Re < 1). The phylodynamic analyses provide epidemiological insights not directly available using standard surveillance and may inform the adjustment of public health strategies in HIV prevention in Ethiopia.

 effective reproductive number effective population size birth–death model phylodynamic HIV-1 epidemic transmission cluster behavioral intervention Ethiopia 2019-05235 2020-02344 Swedish Research Council10.13039/501100004359 Österlund Foundation Medical Faculty, Lund University Max Planck Society10.13039/501100004189 Introduction

The human immunodeficiency virus type 1 (HIV-1) is one of the most devastating infectious diseases in human history (UNAIDS, 2020). At the end of 2020, an estimated 38 million people were living with HIV/AIDS worldwide. Sub-Saharan Africa, the region where HIV-1 emerged during the 1920s, remains the most affected region, accounting for close to 70% of people living with HIV worldwide (Faria et al., 2014; UNAIDS, 2020). Despite the large-scale roll-out of antiretroviral treatment (ART), HIV incidence remains high, mainly in sub-Saharan Africa (UNAIDS, 2020). Ethiopia is one of the many sub-Saharan countries that was severely affected by the HIV epidemic.

HIV-1 is classified into four phylogenetically distinct groups: M (main), N (non-M, non-O), O (outlier), and P (pending), each representing different zoonotic cross-species transmissions of simian immunodeficiency viruses from non-human primates to humans (Sharp and Hahn, 2011; Faria et al., 2014; Giovanetti et al., 2020). Group M is the most prevalent, accounts for more than 95% of all the HIV-1 infections and is divided into 10 subtypes (A–D, F–H, and J–L), more than 102 different circulating recombinant forms (CRFs), and numerous unique recombinant forms (URFs; Hemelaar et al., 2019; Giovanetti et al., 2020). Subtype C is currently the dominant HIV-1 subtype and is responsible for nearly half all HIV-1 infections globally (Hemelaar et al., 2019). Although found worldwide, no official assignment of subtype C strains into phylogenetic sub-subtypes has been made. However, several distinct genetic clades associated with geography have been defined, the southern African clades (C-SA) and the eastern African clade (C-EA). Strains of the C-EA clade and a sub-clade of C-SA, termed C′-ET, are most prevalent in Ethiopia (Thomson and Fernandez-Garcia, 2011; Arimide et al., 2018).

The first HIV-1 infection and AIDS case report in Ethiopia was in 1984 and 1986, respectively (Lester et al., 1988; Tsega et al., 1988). Initially, the epidemic was concentrated to urban areas and along major commercial routes. Serology surveys revealed high prevalence (17%–55%) among risk populations (e.g., female sex workers: FSWs, long-distance truck drivers: LDTD, and soldiers; Mehret, 1990; Mebret et al., 1990). However, after introduction of antiretroviral therapy (ART) in public health care in 2005, the prevalence among the general population decreased and stabilized at significantly lower levels while the prevalence remained high in risk populations (EPHI, 2014).

The HIV epidemic in Ethiopia is considered a generalized epidemic with heterosexual transmission being the dominant mode of transmission (Kebede et al., 2000). Since 1985, Ethiopia has implemented several community-based HIV prevention programs to improve knowledge about the infection and mode of transmission, and interventions to reduce engagement in risk behavior (Mebret et al., 1990; Okubagzhi and Singh, 2002). However, the epidemiological dynamics and their correlations with introduction of various HIV prevention and interventions programs have not been characterized.

Similar to many low-income countries, epidemiological data regarding HIV from Ethiopia are sparse and incomplete, making surveillance of the HIV epidemic challenging. The increased availability of HIV genetic sequencing data and the development of phylogenetic and phylodynamic tools has enabled the use of molecular epidemiology analysis to describe the transmission dynamics and evolutionary history of HIV (Yusim et al., 2001; Delatorre and Bello, 2012; Mir et al., 2018; Vasylyeva et al., 2019). Previous studies in Ethiopia have shown that subtype C dominates the Ethiopian HIV epidemic and have provided valuable insight into HIV genetic diversity, its origins, and epidemic dynamics, but are limited in study participant numbers and geographic and temporal representation (Abebe et al., 2001a,b; Pollakis et al., 2003; Tully and Wood, 2010; Delatorre and Bello, 2012; Mir et al., 2018). Here, we used HIV-1 subtype C pol gene sequences collected from different regions of Ethiopia between 1986 and 2017. We employed state-of-the-art phylogenetic and phylodynamic methods, including both Bayesian coalescent and birth–death modeling, to elucidate evolutionary trajectories and temporal dynamics of the HIV-1 epidemic in Ethiopia.

 Materials and Methods Baseline HIV-1 Drug Resistance Survey

We conducted a prospective HIVDR survey among antiretroviral-naïve adults in St. Paul General Specialized Hospital located in Addis Ababa, Ethiopia, in 2011. We performed the study according to the WHO-recommended survey methodology (Jordan et al., 2008). Treatment-naïve adults (>18 years) eligible to start ART at the St. Paul Generalized Specialized Hospital were consecutively enrolled. Whole blood specimens were collected and transported to the Ethiopia Public Health Institute (EPHI), the national HIV laboratory, and WHO-accredited laboratory for viral load testing and HIVDR genotyping.

HIV genotyping was done using an in-house assay as described previously (Arimide et al., 2018). Briefly, a 1,084 base-pair fragment of HIV-1 pol (corresponding to positions 2,243–3,326 of HXB2; GenBank Accession Number: K03455) comprising amino acids 6–99 of the protease and 1–251 of the reverse transcriptase was obtained by RT-PCR and nested PCR. The purified PCR fragments were then sequenced and analyzed on the ABI 3500xl Genetic Analyzer (Applied Biosystems, Foster City, CA, United States). Sequence assembly and editing were performed using the RECall V 2.0 HIV-1 sequencing analysis tool (University of British Columbia, Vancouver, Canada; Woods et al., 2012). All sequences reported in this study have been deposited in GenBank under Accession Numbers OL598713-OL598856.

 Study Population and Sequence Dataset

We used the dataset of newly sequenced HIV-1 pol sequences and retrieved all publicly available Ethiopian HIV-1 subtype C pol sequences (matching pos. 2,243–3,326 relative of HXB2) from the Los Alamos National Laboratory (LANL) HIV Sequence database1 (Date of access, December 2019). The quality of HIV-1 sequences was verified using the online Quality Control program of the LANL HIV sequence database (see Footnote 1) and sequences with stop codons, frameshifts, and poor quality were removed. We retained only one sequence per patient and selected the earliest sequence for patients with multiple sequences.

We removed duplicate sequences and sequences with potential contamination using the ElimDupes online tool from LANL. Moreover, to identify Ethiopian country-specific transmission clusters, we included a dataset of similar sequences from GenBank by identifying the 10 genetically closest GenBank sequences with BLAST for each Ethiopian HIV-1 subtype C sequence (Altschul et al., 1990; Mount, 2007). We only included sequences of 950 nucleotides or longer with known isolation dates and country of isolation in the analysis, since this 950-bp region has sufficient signal to reconstruct transmission links among infected individuals (Hué et al., 2004).

 HIV-1 Subtyping

Initial explorative HIV-1 subtyping was performed using the online automated subtyping tools REGA v3.0 (Pineda-Peña et al., 2013), COMET v2.2 (Struck et al., 2014), and RIP (Martin et al., 2010). Putative intra-subtype recombinant sequences were detected using jpHMM (jumping profile Hidden Markov Model)2 (Schultz et al., 2009; Arimide et al., 2018). Only non-recombinant sequences were used for the analysis. Final subtyping was determined by maximum likelihood (ML) phylogenetic tree analysis with subtype reference sequences (Arimide et al., 2018).

 Maximum Likelihood Phylogenetic Analyses

A multiple sequence alignment was obtained using MAFFT V. 7 (Katoh and Standley, 2013) and was then manually edited using BioEdit V7.0.9.0 (Hall, 1999) until a non-redundant codon alignment was obtained. To avoid the effect of drug-induced convergent evolution, positions of identified mutations causing or contributing to HIVDR were removed from the alignment, resulting in a final alignment of 909 bp (Wensing et al., 2016).

The initial ML phylogenetic tree was constructed using an online version of PhyML (Guindon et al., 2010) under the GTR + I + Γ4 (general time-reversible nucleotide substitution model using the estimated proportion of invariable sites and four gamma categories). Heuristic tree search was performed using the SPR branch-swapping algorithm. Branch support was determined with aLRT-SH (approximate likelihood ratio test Shimodaira–Hasegawa-like) implemented in PhyML (Guindon et al., 2010). A branch in the phylogeny with an aLRT-SH value ≥0.9 was considered significant (Guindon et al., 2010; Esbjörnsson et al., 2016). The ML trees were visualized using FigTree v1.4.3 (Rambaut, 2016).

Our initial ML phylogenetic trees were constructed using the combined dataset of all Ethiopian sequences and sequences from the BLAST search. To comprehensively describe the HIV-1 subtype C circulating in Ethiopia, the dataset was divided into two based on phylogenetic branch support, the C-EA and C-ET clades.

 Analysis of Transmission Clusters

Separate transmission cluster analysis was performed for the two data sets using the ML phylogenetic analysis implemented by IQ-TREE under GTR + I + Γ4 as selected as the best fitting substitution model for the dataset using jModelTest v2.1.7 and with 1,000 replicates for the aLRT-SH test (Nguyen et al., 2015). Clusters with an aLRT-SH support ≥0.9 were considered significant (Guindon et al., 2010; Esbjörnsson et al., 2016). A transmission cluster was defined as a cluster in the ML phylogeny from root to tips (Esbjörnsson et al., 2016; Hassan et al., 2017; Sallam et al., 2017; Arimide et al., 2018). Clusters with an aLRT-SH-support of ≥0.9 that had a majority (at least 80%) of Ethiopia sequences were considered an Ethiopian transmission cluster. Transmission clusters were also defined based on their sizes (number of sequences/cluster), into dyads (two sequences), medium-sized clusters/networks (3–14 sequences), and large clusters (≥15 sequences; Aldous et al., 2012; Esbjörnsson et al., 2016).

To determine whether there was phylogenetic clustering by geographic region, viral sequences were grouped into six geographic regions (sequence collection location). The strength of association between the geographic location and the phylogeny was determined using two phylogeny–trait association statistics, the parsimony score (PS) and the association index (AI) tests, both of which were implemented in the Bayesian Tip-association Significance testing (BaTS) program (Parker et al., 2008). A significance level of p < 0.05 was used in both statistics.

 Estimating Temporal Signal

For each cluster, we assessed the temporal signal of the data sets by performing root-to-tip genetic distance using TempEst (Rambaut et al., 2016). Clusters that had a positive correlation between genetic diversity and time were considered for further analysis.

 Estimating Viral Phylodynamic History

The birth–death skyline model (BDSKY; Stadler et al., 2012, 2013) implemented in BEAST2 v 2.6.2 was used to quantify epidemic growth through time described by changes in the effective reproductive number (Re) which is the average number of secondary infections from an infected individual at any given time during the epidemic (Bouckaert et al., 2019; Vasylyeva et al., 2019). We used a lognormal distribution prior, LogNorm (0,1), for the effective reproductive number with the upper bound of 10, and a LogNorm (0,1) prior for the become uninfectious rate (δ) in units per year (i.e., the inverse of the time duration of being infectious in a unit of years). We used δ = 0.2, corresponding to a 5-year duration of the infectious period, as the mean of the distribution. In order to account for the uneven number of sequences per year, we employed a different sampling probability (rho) prior for each year, using a beta distribution with a mean equal to the number of samples divided by the reported number of HIV cases in the country for that year. We estimated the change in Re for six equally spaced intervals between the time to most recent common ancestor (tMRCA) and the most recent sampling year.

Phylodynamic analyses were also performed using the Bayesian Skygrid coalescent tree prior, implemented in BEAST 1.10,4 (Gill et al., 2013; Suchard et al., 2018; Hill and Baele, 2019), to estimate changes in effective population size (Ne) through time and estimate the population growth rates (r, years−1) by using a logistic growth coalescent tree prior. Analyses were performed using the GTR + I + Γ4 nucleotide substitution model. The temporal scale of the evolutionary process was estimated using a relaxed uncorrelated molecular clock model with an underlying lognormal distribution with normal priors. This allowed the estimation of the evolutionary rate (μ, nucleotide substitutions per site per year, s/s/y), the age of the most recent common ancestor (tMRCA, years), and the phylodynamic parameters.

For each of the two phylodynamic approaches, we ran three independent Markov Chain Monte Carlo (MCMC) chains until all associated parameters converged to ensure good mixing (ESS > 200) after discarding the first 10% of the MCMC chains. The convergence of the MCMC was inspected visually and by calculating the ESS for each parameter using Tracer v 1.7.5 (Rambaut et al., 2018). We used LogCombiner to combine the different independent results (log and corresponding tree file) from the multiple chains (Drummond et al., 2012). We used the bdskytools package3 in R to plot the results of the BDSKY analysis.

 Ethical Approval

We obtained ethical approval from the Research and Ethical Clearance Committee of EPHI and the National Health Research Ethics Review Committee of the Ministry of Science and Technology of Ethiopia. All participants for the baseline HIV drug resistance survey provided written informed consent to participate in the study.

 Results Study Population and Initial Phylogenetic Analysis

We retrieved 1,132 Ethiopian pol sequences from LANL, collected from different parts of the country from 1986 to 2017. Additionally, we included 144 HIV-1 pol sequences from the baseline HIVDR survey. The combined dataset (n = 1,276 sequences) contained 399 putative recombinant sequences, which were removed from further analysis. We further included a dataset of similar sequences from GenBank by identifying the 10 genetically closest GenBank sequences with BLAST for each of the 877 non-recombinant Ethiopian HIV-1 subtype C sequences in the study. The final combined dataset contained 1,333 non-recombinant HIV-1 subtype C pol sequences (877 Ethiopian and 456 global), which were used for phylogenetic analysis. The ML phylogenetic tree identified two distinct and well-supported clades, the C-EA and C′-ET clades (Figure 1). Among the 877 Ethiopian sequences included in the analysis, the C-EA clade represented 567 (65.0%) of the sequences, while 310 (35.0%) belonged to the C′-ET clade.

Maximum likelihood (ML) phylogenetic tree of HIV-1 subtype C pol sequences (n = 1,333). Maximum likelihood phylogenetic tree constructed using 877 Ethiopian subtype C pol sequences collected 1986–2017 and 456 global subtype C pol sequences. Colored tips are according to the geographic origin of sequences, as indicated in the legend in the top right corner. Branches defining the major clades (C-EA and C′-ET) are those indicated with a filled circle, having branch support (aLRT-SH) >0.9. The Ethiopian clusters are indicated by open circles, having branch support (aLRT-SH) >0.9. The scale bar represents 0.02 substitutions/site.

Most of the Ethiopian C-EA sequences were found in one large cluster (aLRT = 0.87), and only 33 Ethiopian C-EA sequences fell outside this cluster. Sequences of the global dataset intermixed with the Ethiopian sequences and represented sequences obtained most frequently (N = 67, 33.3%) in other East African countries, North America, and Europe. Sequences from Burundi dominated the basally located sequences. In the case of the second major clade, the majority of the Ethiopian C′-ET sequences (95.5%) formed a well-supported sub-clade (aLRT = 0.92), branching off from the basally located sequences. Southern African countries’ sequences were intermixed (N = 121, 47.1%) with the Ethiopian C′-ET sequences, but they were most prominent at the base of the clade.

 Transmission Cluster Analysis

We inferred transmission clusters by separate ML phylogenetic tree analyses of the two clades (C-EA and C′-ET). For C-EA, the ML phylogenetic tree contained a total of 810 sequences (with reference sequences) and identified two large well-supported clusters of 259 (C-EA-259) and 148 sequences (C-EA-148). The C-EA-259 cluster (aLRT = 0.93) contained 213 Ethiopian (82.2% of the sequences of the cluster) and 46 non-Ethiopian sequences, collected 1988–2017, and the C-EA-148 cluster (aLRT = 0.95) contained 124 Ethiopian (83.8%) sequences and 24 non-Ethiopian sequences, collected 1996–2017 (Figure 2). Moreover, we identified six networks (medium-sized clusters) and nine dyads.

Maximum likelihood phylogenetic tree of HIV-1 C-EA clade pol sequences (n = 810). Maximum likelihood phylogenetic tree constructed using 567 Ethiopian subtype C pol sequences collected between 1986 and 2017, 201 global subtype C pol sequences, and 42 reference sequences. Colored tips are according to the geographic origin of sequences, as indicated in the legend in the top left corner. The C-EA-259 and C-EA-148 clusters are highlighted in green, corresponding to a branch support aLRT-SH >0.9 and >80% Ethiopian sequence. The filled circle defining the C′-ET clade represents an aLRT-SH >0.9. The tree was rooted using the C′-ET reference sequences. The scale bar represents 0.03 nucleotide substitutions per site.

Among the 580 C′-ET sequences, we identified one well-supported (aLRT = 0.96) large cluster containing 153 sequences [C′-ET-153; 124 Ethiopian (81.0% of the sequences in the cluster)], and 29 non-Ethiopian sequences collected 1995–2017 (Figure 3). Since we had no information about the associated risk behavior of the respective individuals, we could not associate cluster formation with risk behavior. However, geographic location was available for analysis but was not associated with cluster formation.

Maximum likelihood phylogenetic tree of HIV-1 C′-ET clade pol sequences (n = 580). Maximum likelihood phylogenetic tree constructed using 310 Ethiopian subtype C pol sequences collected 1986–2017, 227 global subtype C pol sequences, and 43 reference sequences. Colored tips are according to the geographic origin of sequences, as indicated in the legend in the top left corner. The C′-ET-153 cluster, highlighted in a green shade, is defined by branch support (aLRT-SH) of >0.9 and >80% Ethiopian sequences. The branches with filled circles possibly defining the C′-ET clade have branch support (aLRT-SH) >0.9. The tree was rooted using the C-EA reference sequences. The scale bar represents 0.02 nucleotide substitutions per site.

 Evolutionary Rates and Dates of HIV-1 Subtype C in Ethiopia

To comprehensively describe the Ethiopian HIV-1 epidemic, we further performed analyses on the three large clusters (C-EA-259, C-EA-148, and C′-ET-153). Root-to-tip analysis indicated a temporal signal in the three data sets (correlation coefficient of 0.53, 0.43, 0.47 for C-EA-259, C-EA-148, and C′-ET-153, respectively).

First, we estimated the tMRCA of the transmission clusters and the C-EA clade. The clade tMRCA represents the date of the origin of the circulating subtype C clade in the region. In contrast, the estimated tMRCA of the Ethiopian transmission clusters should approximate the introductions and local spread of the viral strains in the country (Dalai et al., 2009; Esbjornsson et al., 2011). Based on the inferred tMRCA, the posterior median estimates of C-EA-259 (1975, 95% HPD: 1970–1979) and C-EA-148 (1976, 95% HPD: 1963–1985) were older than the estimate for the C′-ET transmission cluster (1983, 95% HPD: 1975–1988; Table 1). The median root tMRCAs for C-EA clade was estimated to be 1971 (95% HPD: 1966–1976).

Population dynamics and evolutionary estimates for subtype C cluster in Ethiopia.

Subtype C clade/Cluster
C-EA-259 C-EA-148 C′-ET-153
Sequences from Ethiopia (n) 213 124 124
Range of collection (year) 1988–2017 1996–2017 1995–2017
Mean coefficient of variation 0.24 0.23 0.33
Median evolutionary substitution rate (95% HPD)1 1.76 (1.49–2.00) 1.74 (1.11–2.40) 1.83 (1.33–2.32)
Median year of tMRCA (95% HPD) 1975 (1970–1979) 1976 (1963–1985) 1983 (1975–1988)
Median rate of population growth (95% HPD)2 0.66 (0.51–0.81) 0.61 (0.38–0.86) 0.80 (0.53–1.10)
Median epidemic doubling time (years) (95% HPD)3 1.05 (0.86–1.36) 1.12 (0.81–1.82) 0.86 (0.63–1.31)
Maximum effective reproductive number (Re)4 6.13 (95% HPD, 3.53–10.14) 3.93 (95% HPD, 1.88–7.07) 4.88 (95% HPD, 2.57–8.57)
Basic reproductive number (R0)5 4.30 (95% HPD: 3.55–4.05) 4.05 (95% HPD: 2.90–5.25) 5.00 (95% HPD: 3.65–6.50)
Median become uninfectious rate (95% HPD)6 0.13 (0.06, 0.20) 0.19 (0.08, 0.29) 0.21 (0.09, 0.33)

Median number of substitutions/site/year × 10−3.

Median population growth rate (r) per year, determined in BEAST v1.10.4 using a logistic tree prior.

The time (years) required to double the effective number of infections (λ), calculated as λ = ln(2)/r, where r is the population growth rate.

Re (effective reproductive number) which reflect the average number of secondary infections from an infected individual at any given time during the epidemic.

Basic reproductive number (R0) which reflects the average number of infections generated by an infected individual in a population where all individuals are susceptible to infection calculated by using the formula R0 = rD + 1(56) (where r is the population growth rate and D is the average duration of infectiousness period).

Become uninfectious rate, which reflects the inverse of the time duration of being infectious, in the unit of years.

The median estimated evolutionary rate was in the range 1.76–1.82 × 10−3 substitutions/site/year for the three clusters, with overlapping 95% HPD intervals (Table 1), indicating no significant difference of evolutionary rates among the three clusters.

 Temporal Dynamics of Viral Transmission

Direct estimation of the temporal dynamics of the effective reproductive number, Re, was performed using the BDSKY model. The BDSKY analysis assumed a piecewise constant Re, changing over six equidistance intervals between the tMRCA and the most recent sampling. The Re showed similar dynamics for the three clusters. From the start until the beginning of the 1990s, Re remained consistently high (Re > 1) and dropped below the epidemiological threshold (Re < 1) at the mid-1990s and remained below one till recent years (Figures 4A–C). In all three clusters, we observed the maximum Re values (4-6) during the early period (before 1990) of the epidemic (Table 1).

Population dynamics of the HIV-1 epidemic in Ethiopia using the three major clusters (C-EA-259, C-EA-148, and C’ET-153). (A–C) The temporal dynamics of effective reproductive number (Re) using the Bayesian birth–death model, the median Re are shown by the continuous blue line, and indicated in a pink shade is the 95% highest probability density (HPD) intervals. The gray dashed line indicates the last coalescent event reported by the lineage through time (LTT) analysis. The horizontally dotted line represents the epidemiological threshold (Re = 1). (D–F) The median estimates of the effective population size (Ne) over time using the Bayesian skygrid model. The red line shows the median logarithmic effective population size (Ne) over viral generation time (t), representing effective transmissions, and the gray shade indicates the 95% highest probability density (HPD) intervals. The pink dashed line represents the time of antiretroviral therapy (ART) introduction in Ethiopia.

The posterior median estimates of the become non-infectious rate obtained for each cluster were of 0.13 (95% HPD: 0.06–0.20) for C-EA-259, 0.19 (95% HPD: 0.08–0.29) for C-EA-148, and 0.21 (95% HPD: 0.09–0.33) for C′-ET-153, which translates to an infectious period of ~8 years for C-EA-259 and ~5 years for the other two clusters. Despite the longer infectious period estimated for C-EA-259, the overlapping HPD indicates no significant differences among clusters.

We performed different sensitivity analyses to explore the robustness of our BDSKY estimates. We used different values for the mean of the become non-infectious rate prior (δ) going from 6 months to 10 years (δ: 2, 1, 0.5, 0.2, 0.125, and 0.1). We also performed the sensitivity analysis by estimating the effective reproductive number for six and ten equally spaced intervals between tMRCA and the most recent sample. We obtained similar results for all analyses. Re was consistently >1 for the early period until the early 1990s, followed by a decline to Re < 1 after the mid-1990s.

We further performed phylodynamic analysis using the Bayesian Skygrid model to estimate the temporal characteristics of the HIV-1 epidemic in Ethiopia. We analyzed the three Ethiopian clusters and estimated the change in the effective population size (Ne) through time, representing the change in the total number of infections contributing to new cases. The Bayesian skygrid inference revealed a rapid increase in Ne for all the three clusters from the initial introduction period until shortly before the year 2000, followed by a decline and stabilization in Ne until recent years (Figures 4D–F). We also determined the population growth rate (r), the rate of increase in the effective population size with time, using the logistic growth model of the coalescent parametric model. The median growth rate was 0.66, 0.61, and 0.80 year−1 for clusters C-EA-259, C-EA-148, and C′-ET-153, respectively, with overlapping HPD intervals (Table 1).

We also estimated the mean coalescent-based basic reproductive number (R0) values for each cluster from the logistic growth model using the formula R0 = rD + 1 (Pybus et al., 2001; where r is the population growth rate and D is the average duration of infectiousness period). Assuming an average infectious period of 5 years, R0 was in the range 4.0–5.0 for the three clusters, all with overlapping 95% HPD intervals (Table 1).

 Discussion

In this study, we analyzed a large dataset of HIV-1 pol sequences collected from different regions of Ethiopia during more than 30 years. We used both Bayesian coalescent and birth–death models to characterize the dynamics of the HIV-1 epidemic in the country. Overall, our analysis confirms that strains of two subtype C clades are circulating in Ethiopia, supporting the hypothesis that the HIV-1 epidemic in Ethiopia is the result of at least two independent HIV-1 introductions from eastern and southern African countries (Delatorre and Bello, 2012). Moreover, the phylodynamic analyses revealed that the epidemic dynamics in Ethiopia were characterized by an expanding epidemic growth from the start of the epidemic until the mid-1990s, followed by a sharp decline in HIV-1 transmissions. The decline in Re occurred many years before introducing ART and coincided with early behavioral, preventive interventions, and public health awareness campaigns implemented in Ethiopia.

Re is a proxy for HIV incidence and describes the transmission dynamics; Re > 1 means that the epidemic is growing, Re < 1 shows the epidemic is declining, while Re = 1 shows that the epidemic is stabilizing (Stadler et al., 2013). Our phylodynamic analysis showed that the three clusters followed similar epidemic trends. The BDSKY model indicated epidemic growth (Re > 1) from the 1970s to the early 1990s. The basic reproductive number (R0) and mean initial Re were comparably high for the clusters, indicating an early exponential epidemic growth. Similarly, a high epidemic growth rate was estimated for each cluster (0.61–0.80 year−1) and a steady increase in Ne until the beginning of 2000, highlighting the upward trend of HIV transmissions in Ethiopia during the period.

The exponential epidemic growth observed in our analyses is consistent with retrospective serological data, which showed a massive increase of HIV infections among risk populations in Addis Ababa and cities along the main trading routes in Ethiopia during this early period. An extensive survey on FSWs operating in the main trading routes of Ethiopia in 1988 reported an HIV-1 prevalence between 5.3% and 38.1% (Mehret et al., 1990c). Studies performed in the capital Addis Ababa, 1988–90, showed an increase in prevalence from 25% to 54%, and 13% to 18% among FSWs and LDTDs, respectively (Khodakevich et al., 1990; Mehret et al., 1990a,c; Kebede et al., 2000), and 12%–18% among soldiers 1990–1993. Similarly, an increase in HIV prevalence among pregnant women attending antenatal care clinics (ANC) in Addis Ababa (4.6%–10.5%, 1989–90; Kebede et al., 2000) indicated extensive spread in the population.

During the early years, the rapid epidemic increase was most likely due to lack of awareness of HIV, high mobility among FSWs, high-risk sexual behavior, high STI prevalence among the general population (Desta et al., 1990; Mehret et al., 1990b; Negassa et al., 1990), while no prevention interventions were in place. The increased population movement following considerable urbanization and political instability in the country during this early period might also have contributed to the high HIV prevalence and epidemic spread (Hladik et al., 2006; Esbjornsson et al., 2011).

Although there is a lack of data that can describe the HIV epidemic on a national scale, different studies have shown a decline in new infections since the mid-1990s, corroborating our results (Kebede et al., 2000; Tsegaye et al., 2002; Wolday et al., 2007). The decline in HIV prevalence among young adult women (15–24 years) represents a well-established indicator of epidemic decline. It measures the frequency of relatively recent infections and is less influenced by death (Tsegaye et al., 2002). The HIV prevalence trend among young women (15–24 years) attending ANCs in Addis Ababa between 1995 and 2003 declined significantly from 24.2% to 12.9% (Tsegaye et al., 2002; Wolday et al., 2007). Moreover, there was a sharp decline in HIV prevalence among young blood donors in Addis Ababa and nine other towns during this period (Kebede et al., 2000).

Due to a lack of comprehensive data, it has not been easy to obtain estimates of the national incidence trend in Ethiopia. However, a study done to assess the temporal trend among pregnant women who attended the ANCs in the capital Addis Ababa, assessing >7,000 serum specimens collected 1995–2003, showed a significant decline in the HIV-1 incidence rate (from 7.7% to 2.0%, 1995–2003; Wolday et al., 2007). The reduction was substantial among young ANC attendees (aged 15–19 years), indicating an epidemic decline (7.8% to 0.0%, 1995–2003; Wolday et al., 2007). A mathematical modeling study also demonstrated a substantial reduction in the HIV incidence in Ethiopia after 1995 with an estimated annual decline of 6.3% per year, resulting in a total decrease of 77% between 1990 and 2016 (Deribew et al., 2019).

The early decline in the HIV transmissions observed in our study and documented in serological surveys coincide with the change of sexual behavior, prevention, and better control of other sexually transmitted infections (STIs) achieved through the sustained public education and mobilization campaigns. Ethiopia was one of the first countries in sub-Saharan Africa to introduce a task force to prevent and control HIV/AIDS and STI infections, including a national plan for the HIV epidemic response intervention (Zewdie et al., 1990; Kebede et al., 2000; Kloos and Mariam, 2000; Okubagzhi and Singh, 2002). During the early 1990s, Ethiopia had implemented a wide range of HIV prevention and information programs. Implementation of several behavioral interventions and awareness programs took place using the national media, schools, and public gatherings (Hadgu et al., 1990; Zewdie et al., 1990). These programs mainly focused on sustained health education, risk reduction, condom promotion, and prevention and control of STIs (Zewdie et al., 1990; Okubagzhi and Singh, 2002).

The national survey data on behavioral risk factors in Ethiopia are limited. However, different program reviews (1989–1991) and two nationwide surveys on condom use (1987–1993) revealed that these interventions led to changes in sexual risk behavior and increased knowledge about HIV/AIDS. Moreover, the intervention increased condom use and substantially reduced both non-regular partner and STI (Mehret et al., 1996). Similarly, another study showed condom use increased, and non-regular partners decreased among high school students in Addis Ababa and Gondar in the period after 1990 (Kebede et al., 2000). Moreover, a study among male factory workers in Ethiopia showed a change in sexual risk behavior (Mekonnen et al., 2003). Although it is difficult to quantify the impact of the different interventions on HIV incidence, it is reasonable to assume that the various prevention programs impacted HIV transmissions.

Several other studies outside Ethiopia have reported a significant decline in HIV prevalence after behavioral interventions (Martin, 1987; Hessol et al., 1989; Nelson et al., 1996; Stoneburner and Low-Beer, 2004; Halperin et al., 2011). A study in Uganda and Zimbabwe showed a significant decline in HIV prevalence after 1990, resulting from public health intervention on reduced sexual risk behavior (Stoneburner and Low-Beer, 2004; Halperin et al., 2011). Similarly, behavioral interventions resulted in a substantial decrease in HIV transmissions among MSM in Europe and North America in the mid-1980s and heterosexuals in Thailand in the early 1990s (Martin, 1987; Hessol et al., 1989; Nelson et al., 1996; Hué et al., 2005). In line with our results, a comprehensive review of empirical and modeled HIV incidence trends across 20 countries in Sub-Saharan Africa, 1990–2012, revealed a decline in incidence commenced before introducing ART programs, highlighting the significance of behavioral intervention in reducing HIV transmissions (Taaffe et al., 2014).

The trends of the phylodynamic analyses (Figure 4) are in concordance with the UNAIDS HIV incidence and prevalence modeled estimates (Figure 5),4 showing a high incidence and prevalence during the years before 1990–1995, followed by a decline in incidence and stabilization in prevalence.

Modeled mean annual HIV (A) incidence and (B) prevalence, 1990–2018. Data obtained from UNAIDS.

Thus, our results align well with published serological and epidemiological trends in Ethiopia. The epidemic decline coincides with the timing of behavioral interventions in Ethiopia, suggesting a link between the early decline of HIV spread and behavioral interventions many years before the implementation of ART in the country. However, the introduction of ART, which has proved to successfully suppress HIV replication and reduce the risk of onward transmissions, has significantly contributed to reducing HIV transmission, mortality, and maintaining the epidemic decline (Cohen et al., 2011).

The phylogenetic analysis confirms that strains of two HIV-1 subtype C clades (C′-ET and C-EA) are circulating in Ethiopia, suggesting that the HIV epidemic in Ethiopia arose by at least two independent introductions of founder strains from the eastern and southern African countries, respectively (Pollakis et al., 2003; Delatorre and Bello, 2012). Previous studies have defined several distinct subtype C clades that, in most cases, are associated with geographical regions (Thomson and Fernandez-Garcia, 2011). In the case of the Ethiopian lineages, they represent a southern African clade (where the Ethiopian C-SA sub-clade named C′-ET is more or less confined to Ethiopia) and an eastern African clade (C-EA). Our findings align with previous studies showing a distinct phylogeographic subdivision of the HIV-1 subtype C circulating in east, central, and southern African countries (Abebe et al., 2000; Pollakis et al., 2003; Thomson and Fernandez-Garcia, 2011; Delatorre and Bello, 2012).

The transmission cluster analysis indicated that the Ethiopian sequences formed large clusters, indicative of a few major introductions or expansions in the country. The three Ethiopian clusters described here were mixed regarding collection sites, suggesting intermixing of the HIV epidemic in Ethiopia. Different socio-cultural and behavioral factors might also contribute to cluster formation, and assessing these factors are essential for designing HIV-1 transmission preventive strategies. However, the sequences obtained from the public database do not contain information on risk factors, sociodemographic, and other clinical information. Hence, further analysis on factors associated with cluster formation was not possible to discern in this study.

Sequences of Burundi dominated the basal root of the large monophyletic C-EA clade incorporating more than 90% of the Ethiopian sequences, which is in line with a previous study showing that the C-EA clade likely had its origin in Burundi (Delatorre and Bello, 2012). Moreover, the basal root of the monophyletic clade defining the C′-ET clade was dominated by sequences from southern African countries, possibly reflecting the origin of this clade from southern African countries. However, our analysis could not identify the exact countries. Interconnectivity between populations due to geographic proximity has been an essential factor for the spread of HIV across African countries (Wilkinson et al., 2015, 2016; Faria et al., 2019). However, the large distances and cultural interconnectivity between Ethiopia, Burundi, and southern African countries suggest that other factors were in play. Population movements (due to unknown reasons) could have played a role in the introduction of HIV-1 subtype C to Ethiopia, similar to those observed in other parts of Africa (Gray et al., 2009).

Estimating the date of origin and timing of transmissions of HIV is essential to understanding the dynamics of HIV spread. Here, integral to our analysis of HIV transmission dynamics, we also obtained the tMRCA of the C-EA and C′-ET clades in Ethiopia. The molecular dating analysis suggested that the introduction of the C-EA clade took place more than a decade before the first reported AIDS case in Ethiopia. The dating is plausible considering that AIDS symptoms typically arise 6–10 years after infection. Moreover, our tMRCA estimates coincide with estimates of the introductions of the C-EA clade in other Eastern Africa countries, including Kenya, Tanzania, and Uganda (Delatorre and Bello, 2012), and are consistent with previous estimates for subtype C introduction in Ethiopia (Delatorre and Bello, 2012; Mir et al., 2018). Notably, this period also coincided with a large population migration from Burundi, which could have played a crucial role in disseminating the C-EA clade to Eastern Africa countries (Delatorre and Bello, 2012).

In contrast, the tMRCA of C′-ET was estimated at the beginning of the 1980s and is likely the result of a single introduction. This period coincided with the years of socio-political changes in the southern African countries and is associated with a steep growth of the HIV epidemic and viral migrations within southern African countries (Wilkinson et al., 2015, 2016).

To our knowledge, this study represents the most comprehensive study concerning the HIV epidemic in Ethiopia to date. It employs a large number of HIV-1 pol sequences collected during more than 30 years (1986–2017) from different geographical locations in Ethiopia. Moreover, we used state-of-art phylogenetic and phylodynamic methods to investigate the dynamics of the epidemic. Like many other molecular epidemiology studies, we incorporated HIV-1 pol gene sequences deposited in public databases in our analysis. As new HIV infections are recorded, more sequencing will allow to keep track of the ongoing transmission dynamics. The total sampling density was low, mainly due to a generally low sequencing coverage in Ethiopia, compared to the country’s total number of infected individuals. Thus, the transmission clusters identified here cannot fully represent Ethiopia’s entire HIV-1 transmission networks. Moreover, we based our analysis on HIV-1 pol sequences, representing the most sequenced HIV region due to the numerous published HIVDR studies. Although the HIV-1 pol fragment has sufficient phylogenetic signal for phylogenetic analysis of HIV (Hué et al., 2004), longer sequences, including whole genome sequences, may have provided a more informative inference of the HIV-1 molecular epidemiology and transmission history. Finally, the sequences used in this analysis lacked associated information, such as clinical, demographic, risk population assignment, or socio-economic data and, hence, we could not perform a detailed analysis of associated risk factors for HIV transmissions in our study.

In summary, we have employed state-of-art phylogenetic and phylodynamic approaches to describe the molecular epidemiology of HIV in Ethiopia. Our findings indicate that two distinct HIV subtype C strains were introduced in Ethiopia at the beginning of the 1970s and 1980s, followed by rapid epidemic growth until it started to decline in the mid-1990s, a decade before ART roll-out in Ethiopia. The sharp decline coincided with several behavioral prevention interventions and awareness campaigns. Our finding highlights the significance of scaling up behavioral and risk reduction interventions in addition to ART scale-up in the HIV/AIDS control strategy.

 Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

 Author Contributions

DA, PB, and PM conceived and designed the study. DA, YK, PM, and TB coordinated the laboratory tests. DA, LE-G, SS, DK, and PM conducted the phylogenetic and phylodynamic analysis and interpreted the results. DA and PM wrote the manuscript. All authors reviewed the draft and contributed important intellectual content to the final version. All authors agreed and approved to the published version of the manuscript.

 Author Disclaimer

The findings and conclusions in this paper are those of the authors and do not necessarily represent the official position of the funding agencies.

 Funding

PM was supported by the Swedish Research Council (grant numbers 2019-05235 and 2020-02344), Österlund Foundation, and a donation to PM and PB through the Medical Faculty, Lund University, supported the research. LRE-G and DK were supported by the Max Planck Society. The computations were enabled by resources in project SNIC 2021-5-69 provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.821006/full#supplementary-material

The GenBank accession numbers and the year of sampling of sequences used for this study.

We express our deep appreciation to all patients from whom data were collected and analyzed and all staff involved in the HIV-1 drug resistance survey at the St. Paul General Specialized Hospital in Addis Ababa, Ethiopia. We also acknowledge the support given by the EPHI research laboratory team and St. Paul General Specialized Hospital and antiretroviral clinics staff.

 References Abebe A. Lukashov V. V. Pollakis G. Kliphuis A. Fontanet A. L. Goudsmit J. . (2001a). Timing of the HIV-1 subtype C epidemic in Ethiopia based on early virus strains and subsequent virus diversification. AIDS 15, 1555–1561. doi: 10.1097/00002030-200108170-00013, PMID: 11504988 Abebe A. Lukashov V. V. Rinke De Wit T. F. Fisseha B. Tegbaru B. Kliphuis A. . (2001b). Timing of the introduction into Ethiopia of subcluster C' of HIV type 1 subtype C. AIDS Res. Hum. Retrovir. 17, 657–661. doi: 10.1089/088922201300119770, PMID: 11375063 Abebe A. Pollakis G. Fontanet A. L. Fisseha B. Tegbaru B. Kliphuis A. . (2000). Identification of a genetic subcluster of HIV type 1 subtype C (C') widespread in Ethiopia. AIDS Res. Hum. Retrovir. 16, 1909–1914. doi: 10.1089/08892220050195865, PMID: 11118076 Aldous J. L. Pond S. K. Poon A. Jain S. Qin H. Kahn J. S. . (2012). Characterizing HIV transmission networks across the United States. Clin. Infect. Dis. 55, 1135–1143. doi: 10.1093/cid/cis612, PMID: 22784872 Altschul S. F. Gish W. Miller W. Myers E. W. Lipman D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2 Arimide D. A. Abebe A. Kebede Y. Adugna F. Tilahun T. Kassa D. . (2018). HIV-genetic diversity and drug resistance transmission clusters in Gondar, northern Ethiopia, 2003-2013. PLoS One 13:e0205446. doi: 10.1371/journal.pone.0205446, PMID: 30304061 Bouckaert R. Vaughan T. G. Barido-Sottani J. Duchêne S. Fourment M. Gavryushkina A. . (2019). BEAST 2.5: an advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15:e1006650. doi: 10.1371/journal.pcbi.1006650, PMID: 30958812 Cohen M. S. Chen Y. Q. McCauley M. Gamble T. Hosseinipour M. C. Kumarasamy N. . (2011). Prevention of HIV-1 infection with early antiretroviral therapy. N. Engl. J. Med. 365, 493–505. doi: 10.1056/NEJMoa1105243, PMID: 21767103 Dalai S. C. de Oliveira T. Harkins G. W. Kassaye S. G. Lint J. Manasa J. . (2009). Evolution and molecular epidemiology of subtype C HIV-1 in Zimbabwe. AIDS 23, 2523–2532. doi: 10.1097/QAD.0b013e3283320ef3, PMID: 19770693 Delatorre E. O. Bello G. (2012). Phylodynamics of HIV-1 subtype C epidemic in East Africa. PLoS One 7:e41904. doi: 10.1371/journal.pone.0041904, PMID: 22848653 Deribew A. Biadgilign S. Deribe K. Dejene T. Tessema G. A. Melaku Y. A. . (2019). The burden of HIV/AIDS in Ethiopia from 1990 to 2016: evidence from the global burden of diseases 2016 study. Ethiop. J. Health Sci. 29, 859–868. doi: 10.4314/ejhs.v29i1.7, PMID: 30700953 Desta S. Feleke W. Yusuf M. Mehiret M. Geyid A. Ghidinelli M. . (1990). Prevalence of STD and STD related risk factors in sex workers of Addis Ababa. Ethiop. J. Heal. Dev. 4, 149–153. Drummond A. J. Suchard M. A. Xie D. Rambaut A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. doi: 10.1093/molbev/mss075, PMID: 22367748 EPHI (2014). Ethiopian National Key Population HIV Bio-Behavioral Surveillance Round I, 2013 Report: EPHI; 2014. Esbjörnsson J. Mild M. Audelin A. Fonager J. Skar H. Bruun Jørgensen L. . (2016). HIV-1 transmission between MSM and heterosexuals, and increasing proportions of circulating recombinant forms in the Nordic countries. Virus Evol. 2:vew010. doi: 10.1093/ve/vew010, PMID: 27774303 Esbjornsson J. Mild M. Mansson F. Norrgren H. Medstrand P. (2011). HIV-1 molecular epidemiology in Guinea-Bissau, West Africa: origin, demography and migrations. PLoS One 6:e17025. doi: 10.1371/journal.pone.0017025, PMID: 21365013 Faria N. R. Rambaut A. Suchard M. A. Baele G. Bedford T. Ward M. J. . (2014). HIV epidemiology. The early spread and epidemic ignition of HIV-1 in human populations. Science 346, 56–61. doi: 10.1126/science.1256739, PMID: 25278604 Faria N. R. Vidal N. Lourenco J. Raghwani J. Sigaloff K. C. E. Tatem A. J. . (2019). Distinct rates and patterns of spread of the major HIV-1 subtypes in central and East Africa. PLoS Pathog. 15:e1007976. doi: 10.1371/journal.ppat.1007976, PMID: 31809523 Gill M. S. Lemey P. Faria N. R. Rambaut A. Shapiro B. Suchard M. A. (2013). Improving Bayesian population dynamics inference: a coalescent-based model for multiple loci. Mol. Biol. Evol. 30, 713–724. doi: 10.1093/molbev/mss265, PMID: 23180580 Giovanetti M. Ciccozzi M. Parolin C. Borsetti A. (2020). Molecular epidemiology of HIV-1 in African countries: a comprehensive overview. Pathogens 9:1072. doi: 10.3390/pathogens9121072, PMID: 33371264 Gray R. R. Tatem A. J. Lamers S. Hou W. Laeyendecker O. Serwadda D. . (2009). Spatial phylodynamics of HIV-1 epidemic emergence in East Africa. AIDS 23, F9–F17. doi: 10.1097/QAD.0b013e32832faf61, PMID: 19644346 Guindon S. Dufayard J. F. Lefort V. Anisimova M. Hordijk W. Gascuel O. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59, 307–321. doi: 10.1093/sysbio/syq010, PMID: 20525638 Hadgu T. G. Egziabher E. Gizaw G. Yilma A. Khodakevich L. Zewdie D. . (1990). Intersectoral collaboration in AIDS control in Ethiopia. Ethiop. J. Heal. Dev. 4, 7–99. Hall T. A. (1999). BioEdit: a user-friendly biological sequence alignment editor and analysis program for windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98. Halperin D. T. Mugurungi O. Hallett T. B. Muchini B. Campbell B. Magure T. . (2011). A surprising prevention success: why did the HIV epidemic decline in Zimbabwe? PLoS Med. 8:e1000414. doi: 10.1371/journal.pmed.1000414, PMID: 21346807 Hassan A. S. Pybus O. G. Sanders E. J. Albert J. Esbjörnsson J. (2017). Defining HIV-1 transmission clusters based on sequence data. AIDS 31, 1211–1222. doi: 10.1097/QAD.0000000000001470, PMID: 28353537 Hemelaar J. Elangovan R. Yun J. Dickson-Tetteh L. Fleminger I. Kirtley S. . (2019). Global and regional molecular epidemiology of HIV-1, 1990-2015: a systematic review, global survey, and trend analysis. Lancet Infect. Dis. 19, 143–155. doi: 10.1016/S1473-3099(18)30647-9, PMID: 30509777 Hessol N. A. Lifson A. R. O'Malley P. M. Doll L. S. Jaffe H. W. Rutherford G. W. (1989). Prevalence, incidence, and progression of human immunodeficiency virus infection in homosexual and bisexual men in hepatitis B vaccine trials, 1978-1988. Am. J. Epidemiol. 130, 1167–1175. doi: 10.1093/oxfordjournals.aje.a115445, PMID: 2531543 Hill V. Baele G. (2019). Bayesian estimation of past population dynamics in BEAST 1.10 using the Skygrid coalescent model. Mol. Biol. Evol. 36, 2620–2628. doi: 10.1093/molbev/msz172, PMID: 31364710 Hladik W. Shabbir I. Jelaludin A. Woldu A. Tsehaynesh M. Tadesse W. (2006). HIV/AIDS in Ethiopia: where is the epidemic heading? Sex. Transm. Infect. 82(Suppl. 1), i32–i35. doi: 10.1136/sti.2005.016592, PMID: 16581757 Hué S. Clewley J. P. Cane P. A. Pillay D. (2004). HIV-1 pol gene variation is sufficient for reconstruction of transmissions in the era of antiretroviral therapy. AIDS 18, 719–728. doi: 10.1097/00002030-200403260-00002, PMID: 15075506 Hué S. Pillay D. Clewley J. P. Pybus O. G. (2005). Genetic analysis reveals the complex structure of HIV-1 transmission within defined risk groups. Proc. Natl. Acad. Sci. U. S. A. 102, 4425–4429. doi: 10.1073/pnas.0407534102, PMID: 15767575 Jordan M. R. Bennett D. E. Bertagnolio S. Gilks C. F. Sutherland D. (2008). World Health Organization surveys to monitor HIV drug resistance prevention and associated factors in sentinel antiretroviral treatment sites. Antivir. Ther. 13(Suppl. 2), 15–23. doi: 10.1186/s13104-016-2101-8 Katoh K. Standley D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. doi: 10.1093/molbev/mst010, PMID: 23329690 Kebede D. Aklilu M. Sanders E. (2000). The HIV epidemic and the state of its surveillance in Ethiopia. Ethiop. Med. J. 38, 283–302. PMID: 11125503 Khodakevich L. Mehret M. Negassa H. Shanko B. (1990). Progression of human immunodeficeiencyy virus epidemic in Ethiopia. Ethiop. J. Heal. Dev. 4, 183–187. Kloos H. Mariam D. H. (2000). HIV/AIDS in Ethiopia: an overview. Northeast. Afr. Stud. 7, 13–40. doi: 10.1353/nas.2004.0006 Lester F. T. Ayehunie S. Zewdie D. (1988). Acquired immunodeficiency syndrome: seven cases in an Addis Ababa hospital. Ethiop. Med. J. 26, 139–145. PMID: 3416846 Martin J. L. (1987). The impact of AIDS on gay male sexual behavior patterns in New York City. Am. J. Public Health 77, 578–581. doi: 10.2105/AJPH.77.5.578, PMID: 3565650 Martin D. P. Lemey P. Lott M. Moulton V. Posada D. Lefeuvre P. (2010). RDP3: a flexible and fast computer program for analyzing recombination. Bioinformatics 26, 2462–2463. doi: 10.1093/bioinformatics/btq467, PMID: 20798170 Mebret M. Kbodakcvicb L. Zewdie D. Ayebunic S. Gizaw G. Shanko B. . (1990). HIV-1 infection and related risk factors among female sex workers in urban areas of Ethiopia. Ethiop. J. Heal. Dev. 4, 163–170. Mehret M. (1990). HIV-1 infection and related risk factors among female sex workers in urban areas in Ethiopia. Ethiop. J. Health Dev. 4(Suppl. 2), 163–170. Mehret M. Khodakevich L. Zewdie D. Gizaw G. Ayehune S. Shanko B. . (1990a). HIV-1 infection among employees of the Ethiopian Freightt transport corporation. Ethiop. J. Heal. Dev. 4, 177–182. Mehret M. Mertens T. E. Caraël M. Negassa H. Feleke W. Yitbarek N. . (1996). Baseline for the evaluation of an AIDS programme using prevention indicators: a case study in Ethiopia. Bull. World Health Organ. 74, 509–516. PMID: 9002331 Mehret M. K. L. Shanko B. Belete F. (1990b). Sexual behaviours and some social features off female sex workers in the city of Addis Ababa. Ethiop. J. Health Dev. 4, 133–113. Mehret M. K. L. Zewdie D. Ayehunie S. Shanko B. Gizaw G. . (1990c). HIV-1 infection and some related risk factors among female sex workers in Addis Ababa. Ethiop. J. Health Dev. 4, 171–176. Mekonnen Y. Sanders E. Aklilu M. Tsegaye A. Rinke de Wit T. F. Schaap A. . (2003). Evidence of changes in sexual behaviours among male factory workers in Ethiopia. AIDS 17, 223–231. doi: 10.1097/00002030-200301240-00013, PMID: 12545083 Mir D. Graf T. de Matos E. Almeida S. Pinto A. R. Delatorre E. . (2018). Inferring population dynamics of HIV-1 subtype C epidemics in eastern Africa and southern Brazil applying different Bayesian phylodynamics approaches. Sci. Rep. 8:8778. doi: 10.1038/s41598-018-26824-4, PMID: 29884822 Mount D. W. (2007). Using the basic local alignment search tool (BLAST). Cold Spring Harb. Protoc. 2007:pdb.top17. doi: 10.1101/pdb.top17 Negassa H. Kefene H. Khodakevich L. Zewdie D. Shanko B. (1990). Profile of AIDS case in Ethiopia. Ethiop. J. Heal. Dev. 4, 213–217. Nelson K. E. Celentano D. D. Eiumtrakol S. Hoover D. R. Beyrer C. Suprasert S. . (1996). Changes in sexual behavior and a decline in HIV infection among young men in Thailand. N. Engl. J. Med. 335, 297–303. doi: 10.1056/NEJM199608013350501 Nguyen L. T. Schmidt H. A. von Haeseler A. Minh B. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274. doi: 10.1093/molbev/msu300, PMID: 25371430 Okubagzhi G. Singh S. (2002). Establishing an HIV/AIDS programme in developing countries: the Ethiopian experience. AIDS 16, 1575–1586. doi: 10.1097/00002030-200208160-00002, PMID: 12172078 Parker J. Rambaut A. Pybus O. G. (2008). Correlating viral phenotypes with phylogeny: accounting for phylogenetic uncertainty. Infect. Genet. Evol. 8, 239–246. doi: 10.1016/j.meegid.2007.08.001, PMID: 17921073 Pineda-Peña A. C. Faria N. R. Imbrechts S. Libin P. Abecasis A. B. Deforche K. . (2013). Automated subtyping of HIV-1 genetic sequences for clinical and surveillance purposes: performance evaluation of the new REGA version 3 and seven other tools. Infect. Genet. Evol. 19, 337–348. doi: 10.1016/j.meegid.2013.04.032, PMID: 23660484 Pollakis G. Abebe A. Kliphuis A. De Wit T. F. Fisseha B. Tegbaru B. . (2003). Recombination of HIV type 1C (C′/C″) in Ethiopia: possible link of EthHIV-1C′ to subtype C sequences from the high-prevalence epidemics in India and Southern Africa. AIDS Res. Hum. Retrovir. 19, 999–1008. doi: 10.1089/088922203322588350, PMID: 14678607 Pybus O. G. Charleston M. A. Gupta S. Rambaut A. Holmes E. C. Harvey P. H. (2001). The epidemic behavior of the hepatitis C virus. Science 292, 2323–2325. doi: 10.1126/science.1058321 Rambaut A. (2016). FigTree v1.4.3: Tree Figure Drawing Tool. Available at: http://tree.bio.ed.ac.uk/software/figtree/ (Accessed October 4, 2019). Rambaut A. Drummond A. J. Xie D. Baele G. Suchard M. A. (2018). Posterior summarization in Bayesian Phylogenetics using tracer 1.7. Syst. Biol. 67, 901–904. doi: 10.1093/sysbio/syy032, PMID: 29718447 Rambaut A. Lam T. T. Max Carvalho L. Pybus O. G. (2016). Exploring the temporal structure of heterochronous sequences using TempEst (formerly path-O-gen). Virus Evol. 2:vew007. doi: 10.1093/ve/vew007, PMID: 27774300 Sallam M. Esbjornsson J. Baldvinsdottir G. Indriethason H. Bjornsdottir T. B. Widell A. . (2017). Molecular epidemiology of HIV-1 in Iceland: early introductions, transmission dynamics and recent outbreaks among injection drug users. Infect. Genet. Evol. 49, 157–163. doi: 10.1016/j.meegid.2017.01.004, PMID: 28082188 Schultz A. K. Zhang M. Bulla I. Leitner T. Korber B. Morgenstern B. . (2009). jpHMM: improving the reliability of recombination prediction in HIV-1. Nucleic Acids Res. 37, W647–W651. doi: 10.1093/nar/gkp371, PMID: 19443440 Sharp P. M. Hahn B. H. (2011). Origins of HIV and the AIDS pandemic. Cold Spring Harb. Perspect. Med. 1:a006841. doi: 10.1101/cshperspect.a006841 Stadler T. Kouyos R. von Wyl V. Yerly S. Boni J. Burgisser P. . (2012). Estimating the basic reproductive number from viral sequence data. Mol. Biol. Evol. 29, 347–357. doi: 10.1093/molbev/msr217, PMID: 21890480 Stadler T. Kuhnert D. Bonhoeffer S. Drummond A. J. (2013). Birth-death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). Proc. Natl. Acad. Sci. 110, 228–233. doi: 10.1073/pnas.1207965110, PMID: 23248286 Stoneburner R. L. Low-Beer D. (2004). Population-level HIV declines and behavioral risk avoidance in Uganda. Science 304, 714–718. doi: 10.1126/science.1093166, PMID: 15118157 Struck D. Lawyer G. Ternes A. M. Schmit J. C. Bercoff D. P. (2014). COMET: adaptive context-based modeling for ultrafast HIV-1 subtype identification. Nucleic Acids Res. 42:e144. doi: 10.1093/nar/gku739, PMID: 25120265 Suchard M. A. Lemey P. Baele G. Ayres D. L. Drummond A. J. Rambaut A. (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4:vey016. doi: 10.1093/ve/vey016, PMID: 29942656 Taaffe J. Fraser-Hurt N. Gorgens M. Harimurti P. (2014). A Comprehensive Review of Empirical and Modeled HIV Incidence Trends (1990–2012). Policy Research Working Paper; No. 7042. World Bank. Available at: http://hdl.handle.net/10986/20378 (Accessed March 10, 2021). Thomson M. M. Fernandez-Garcia A. (2011). Phylogenetic structure in African HIV-1 subtype C revealed by selective sequential pruning. Virology 415, 30–38. doi: 10.1016/j.virol.2011.03.021, PMID: 21507449 Tsega E. Mengesha B. Nordenfelt E. Hansson B. G. Lindberg J. (1988). Serological survey of human immunodeficiency virus infection in Ethiopia. Ethiop. Med. J. 26, 179–184. PMID: 3215178 Tsegaye A. Rinke De Wit T. F. Mekonnen Y. Beyene A. Aklilu M. Messele T. . (2002). Decline in prevalence of HIV-1 infection and syphilis among young women attending antenatal care clinics in Addis Ababa, Ethiopia: results from sentinel surveillance, 1995-2001. J. Acquir. Immune Defic. Syndr. 30, 359–362. doi: 10.1097/00126334-200207010-00013, PMID: 12131574 Tully D. C. Wood C. (2010). Chronology and evolution of the HIV-1 subtype C epidemic in Ethiopia. AIDS 24, 1577–1582. doi: 10.1097/QAD.0b013e32833999e1, PMID: 20539092 UNAIDS (2020). Global HIV & AIDS statistics—Fact Sheet. Available at: https://www.unaids.org/en/resources/fact-sheet (Accessed March 10, 2021). Vasylyeva T. I. du Plessis L. Pineda-Pena A. C. Kuhnert D. Lemey P. Vandamme A. M. . (2019). Tracing the impact of public health interventions on HIV-1 transmission in Portugal using molecular epidemiology. J. Infect. Dis. 220, 233–243. doi: 10.1093/infdis/jiz085, PMID: 30805610 Wensing A. M. Calvez V. Günthard H. F. Johnson V. A. Paredes R. Pillay D. . (2016). 2017 update of the drug resistance mutations in HIV-1. Top. Antivir. Med. 24, 132–133. PMID: 28208121 Wilkinson E. Engelbrecht S. de Oliveira T. (2015). History and origin of the HIV-1 subtype C epidemic in South Africa and the greater southern African region. Sci. Rep. 5:16897. doi: 10.1038/srep16897, PMID: 26574165 Wilkinson E. Rasmussen D. Ratmann O. Stadler T. Engelbrecht S. de Oliveira T. (2016). Origin, imports and exports of HIV-1 subtype C in South Africa: a historical perspective. Infect. Genet. Evol. 46, 200–208. doi: 10.1016/j.meegid.2016.07.008 Wolday D. Meles H. Hailu E. Messele T. Mengistu Y. Fekadu M. . (2007). Temporal trends in the incidence of HIV infection in antenatal clinic attendees in Addis Ababa, Ethiopia, 1995-2003. J. Intern. Med. 261, 132–137. doi: 10.1111/j.1365-2796.2006.01740.x, PMID: 17241178 Woods C. K. Brumme C. J. Liu T. F. Chui C. K. Chu A. L. Wynhoven B. . (2012). Automating HIV drug resistance genotyping with RECall, a freely accessible sequence analysis tool. J. Clin. Microbiol. 50, 1936–1942. doi: 10.1128/JCM.06689-11, PMID: 22403431 Yusim K. Peeters M. Pybus O. G. Bhattacharya T. Delaporte E. Mulanga C. . (2001). Using human immunodeficiency virus type 1 sequences to infer historical features of the acquired immune deficiency syndrome epidemic and human immunodeficiency virus evolution. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 356, 855–866. doi: 10.1098/rstb.2001.0859, PMID: 11405933 Zewdie D. Gizaw G. Khodakevich L. Degifie G. Wemeue M. (1990). Development and management of the AIDS control programme in Ethiopia. Ethiop. J. Heal. Dev. 4, 87–96.

1http://www.hiv.lanl.gov

2http://jphmm.gobics.de/submission_hiv

3https://github.com/laduplessis/bdskytools

4https://aidsinfo.unaids.org