Study design and site selection

A cross-sectional survey was conducted between August 2011 and December 2013 among antiretroviral-naive adults to evaluate transmitted drug resistance, according to the World Health Organization (WHO)-recommended threshold survey methodology [13]. It was conducted in two of the major Voluntary Counselling and Testing (VCT) clinics in Gondar, located 700 km north of Addis Ababa. Gondar is the second largest town in Amhara region and, it was one of the areas in the country where public ART was first initiated (in year 2003). The HIV-prevalence 2003–2013 in Gondar was on a stable level (mean: 10.6%) while the national prevalence during the same time period showed a declining trend, from 12% to 4.4% [32]. Furthermore, three studies conducted at the hospital in 2003–2010 showed the presence of HIVDR in the area [23–25]. The sample size for the current survey followed the sequential sampling method selected by WHO for the surveillance of transmitted HIVDR in low-resource settings. According to the recommendation, it is advised to collect about 70 specimens of eligible individuals consecutively diagnosed with HIV in sites within a survey area [20].

Study participants

Individuals with new diagnosed HIV infection, among VCT clients at the two survey sites were asked to participate in the study. By adopting the WHO recommended inclusion criterion, participants 18–25 years old, who had no prior history of HIV/AIDS-related illness and no history of ART, resident of Gondar for more than a year, and had no history of previous pregnancy were consecutively enrolled. After obtaining written informed consent, 10 ml of blood was collected by venepuncture. Screening for HIV was done using point of care rapid testing format employed for HIV diagnosis in Ethiopia. This algorithm uses HIV (1 + 2) Antibody Colloidal Gold (KHB, Shanghai Kehua Bio-engineering Co Ltd, China) as a screening test, followed by HIV 1/2 STAT-PAK (Chembio Diagnostics, USA) if positive. In cases with negative STAT-PAK results following a positive KHB test, a third test, Unigold HIV (Trinity Biotech, Ireland), was used as confirmation. Specimens were transported to Ethiopia Public Health Institute, National HIV laboratory, Addis Ababa (a WHO accredited laboratory for genotyping) on dry ice for long term storage at -80°C until genotyping. At the time of blood sampling basic demographic and clinical information, including age, gender, and history of HIV test was collected using a standardized questionnaire.

HIV-1 genotyping

A 1084-bp fragment of HIV-1 pol (corresponding to the position: 2243–3326 of HXB2, Genbank Accession Number: K03455) comprising amino acids 6–99 of the protease (PR) and 1–251 of the reverse transcriptase (RT) was amplified using an in-house genotyping assay as described in [33], and S1 Text. PCR products were directly sequenced using six primers (three on each strand) on an ABI 3100 or an ABI 3500xl DNA Genetic Analyzer (Applied Biosystems). Sequence assembly and editing were performed using RECall V 2.0 HIV-1 sequencing analysis tool [34].

Identification of drug resistance mutations

Surveillance drug resistance mutations (SDRMs) were examined according to the Stanford Genotypic Resistance calibrated population resistance (CPR) tool version 6.0 based on the WHO surveillance transmitted drug resistance mutation list of 2009 [35, 36]. Classification of TDR level (low: < 5%, moderate: 5–15%, or high: >15%) was made based on the WHO threshold survey protocol [20].

HIV-1 subtyping and recombination analysis

Sequence quality control was performed using the online Quality Control program of the Los Alamos HIV sequence database (hiv.lanl.gov). The REGA and Comet online subtyping tools were used for initial classification into subtypes and inter-subtype recombinants [37, 38]. Putative intra-subtype recombinants were verified using Simplot ver. 3.5.1[39]. All sequences were also screened for recombination using RDP ver. 3.44 [40]. Final subtyping was performed through phylogenetic analysis using the reference HIV-1 data set from Los Alamos HIV sequence database (hiv.lanl.gov). Sequences were aligned using ClustalX2 [41] and then edited to a final length of 1044 bases using BioEdit v4.0.6 (http://www.mbio.ncsu.edu/bioedit/bioedit.html). A Maximum likelihood (ML) phylogenetic tree was constructed using the online version of PhyML with the GTR+I+Γ nucleotide substitution model (using estimated proportion of invariable sites and four gamma categories) and NNI plus SPR to estimate the tree topology. Branch support was determined with aLRT-SH (approximate likelihood ratio test Shimodaira-Hasegawa like) implemented in PhyML. A branch in the phylogeny with an aLRT-SH value ≥0.9 was considered significant [42, 43].

To further dissect the subtype C distribution in Gondar, sub-subtyping into subtype C clades and detection of putative intra-subtype C recombinants were performed using sequences obtained from this study population and all previously reported HIV-1 pol sequences from Gondar [23–25]. Briefly, reference sequences of the three subtype C clades were obtained as described in S1 Text. Putative intra-subtype recombinant subtype C sequences were identified by jpHMM [44] using parameters as outlined in S1 Text. The final subtyping ML tree is shown in S1 Fig.

Clade distribution and transmission cluster analysis

For the purpose of a comprehensive analysis of temporal changes of subtype C clades over time in Gondar, subtype C clades were identified using ML phylogenetic tree analysis. For the purpose of transmission cluster analysis, we followed a methodology described previously and included a data set of similar sequences from GenBank by identifying the ten best scoring GenBank sequences with BLAST for each of the Gondar sequences in the study [43, 45, 46]. The final data set contained 491 sequences (n = 299 from Gondar and n = 192 GenBank reference sequences; see S2 Table) which were used to define local (Gondar) transmission clusters. Transmission clusters were defined as described previously [43, 45–48]. Briefly, a transmission cluster was defined as a cluster in the ML phylogeny from root to tips. A cluster with an aLRT SH-support of ≥0.9 that had a majority (at least 80%) of sequences from Gondar was considered as a Gondar (local) transmission cluster. Transmission clusters were 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) [43, 49]. Clusters sharing >33% of the same DRM were defined as a putative drug resistance transmission clusters [50].

Evolutionary and phylodynamic analysis

To estimate an evolutionary rate for each of the three clades which were needed to perform dated cluster analysis using phylodynamic analyses, three data sets were assembled, containing both global and Ethiopian sequences which were considered to be representative for the different clades. This approach has been used previously for HIV-1 subtype B [51, 52]. The analysis was performed using the sequences shown in S1 Fig and S2 Table, by randomly selecting a subset of sequences from each clade (N = 86, N = 65 and N = 70 for C-EA, C’-ET and C-SA, respectively). A maximum of one sequence of each transmission cluster was allowed in the final data set. For the data set containing the southern African clade (C-SA) sequences, the temporal signal was weak (R² = 0.03), as assessed by root-to-tip analysis using TempEst, indicating that the data set was not optimal for estimating reliable substitution rates, while the R² was higher for the C-EA and C’-ET clades (R² = 0.35 and 0.22, respectively), indicating a better temporal signal in the datasets [53]. We therefore estimated the substitution rates for the C-EA and C’-ET sequences only and relied on previous estimates for the C-SA clades (see below). For these analyses we centered the prior mean rate at 0.001 substitutions/site/year and specified the standard deviation to 0.33 on a lognormal distribution such that the 2.5 and 97.5 percentiles of the distribution contained the rates obtained previously of HIV-1 pol of subtype C [54, 55]. The evolutionary rate was estimated by employing the Bayesian Markov Chains Monte Carlo (MCMC) clock method implemented in the BEAST software package v1.8.4 [56]. We estimated the evolutionary rate using both a strict and relaxed Bayesian MCMC clock, in both cases with a flexible demographic model (the Bayesian skyline plot) as a tree prior. This analysis estimated posterior distribution of median rates for the C-EA and C’-ET clades to 1.26 x 10⁻³ (95% HPD: 8.80x10^-4–1.64x10^-3) and 1.44x10^-3 (95% HPD: 9.70x10^-4–1.91x10^-3) substitutions/site/year, respectively (S1 Table). Using these rates as priors, subsequent BEAST analysis was performed on the C-EA and C’-ET sequences obtained from Gondar only while a previous evolutionary rate of 2.15x10-3 (1.79x10^-3–2.60x10^-3, 95% HPD) substitutions/site/year was used as a prior for the C-SA sequences of Gondar [54]. In all cases, a lognormal distribution on the rate prior was employed, as described above, such that the prior mean rate was centered at each clade-specific mean rate and that the 2.5 and 97.5 percentiles of the distribution contained the rates obtained from the above analyses (S1 Table). For all analyses we used both strict and relaxed clock models and specified a codon position partitioning model (the SRD06 model) for all data sets [57]. We also employed different demographic tree priors on each data set, the skyride, the logistic and exponential demographic models implemented in BEAST (S1 Table). Model comparisons were done by Bayes factor (BF) analysis of marginal likelihoods [58]. BF values larger than 3 and 5 were considered to represent strong and very strong evidence respectively, against H₀ [59].

Statistical analysis

Statistical tests were performed using SPSS 24 (IBM Corp., Armonk, NY, USA). DR prevalence was determined with a confidence interval (CI) of 95% using the Wilson method. Categorical variables were compared using 2-tailed Fisher’s exact test, while continuous variables were compared using Mann-Whitney 2-tailed U test. Trends over time were analysed using linear-by-linear test for association.

Ethical approval

Scientific and ethical approval was granted by the Research and Ethical Clearance Committee of the Ethiopian Public Health Institute, and the National Health Research Ethics Review Committee of Ministry of Science and Technology of Ethiopia. All participants provided written informed consent.

Availability of data

Nucleotide sequences reported in this study have been deposited in the Genbank repository (Accession Numbers: KM390990—KM391398).

Population characteristics

A total of 84 individuals attending the VCT clinics at Gondar University hospital (n = 64) and Gondar health centre (n = 20) met the inclusion criteria and were enrolled in the study.

All samples were confirmed as HIV-1 positive. In total, 78% (67/84) of the samples were successfully amplified and sequenced. The median age among the 84 individuals was 22 years (IQR: 20–23) and 83% were females, which was similar to the age (22 years [IQR: 20–24]) and gender distribution (87% females) of the 67 participants whose samples were successfully genotyped (p = 0.992 and p = 0.653, M-W and Fisher exact test, respectively).

Levels of transmitted drug resistance

TDR levels were determined from the 67 sequenced samples. Four specimens (6.0%) carried major DRMs, corresponding to a moderate TDR level [60]. Two had K103N, one G190S and one Y181C, all representing resistance to non-nucleoside reverse transcriptase inhibitors (NNRTI). No NRTI or protease inhibitor (PI) associated DRMs were identified among the sequenced samples (Table 1). Although there were more females enrolled in the survey than men there were no association between gender and DRMs (2 men and 2 females carried HIVDR virus; p = 0.084, FET).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Prevalence of HIV drug resistance mutations in Gondar 2003–2013.

https://doi.org/10.1371/journal.pone.0205446.t001

Level of HIV-1 drug resistance in Gondar 2003–2013

Analysis of HIV DRMs of two previous studies from Gondar after ART roll out in 2009 and 2010 as described in Methods, revealed levels comparable to our threshold study (4.4–8.2% versus 6.0%; Table 1) [23, 25]. These estimates were higher than the levels reported in the study performed before ART roll-out in Gondar 2003 (3.3%) [24]. However, the confidence intervals of the levels of DRMs were largely overlapping between the studies and the proportion of DRMs before and after ART roll out were not different (p = 0.584, FET; Table 1). DRMs of the previous three studies were as in the present study, most often associated with NNRTIs. Taken together, the G190A/S/E mutations were most prominent (seven of the 12 NNRTI mutations), followed by the K103N, Y181I/C and K101E substitutions (two, two and one, respectively). Mutations conferring resistance to NRTIs and PIs were found in earlier studies at low levels (1.0–1.9%) but were not identified in the present study (Table 1).

Subtype C clade distribution and temporal changes 2003–2013

Phylogenetic subtyping revealed that 94% (63 of 67) of study participants enrolled 2011–2013 were infected with HIV-1 subtype C, while the remaining four study participants had subtype A, B or A/C recombinant sequences (S2 Table). Since previous observations have suggested more rapid expansion of the C’-ET strains than the C-EA strains in Ethiopia, we further investigated the distribution of our sequences into different subtype C clades [61]. Detailed analysis indicated that 59 of 63 sequences represented the three major subtype C clades: C-EA (n = 30), C’-ET (n = 16) and C-SA (n = 13), while four sequences were putative C-EA/C’-ET recombinants (S2 Table). For in-depth phylogenetic analysis, we next constructed a new data set that represented non-recombinant subtype C sequences of the current study (n = 59), the three previous studies from Gondar in 2003–2010 (n = 242) and a Genbank reference data set (n = 190), resulting in a final data set of 491 sequences collected 1986–2013 (n = 301 from Gondar; S2 Table). Phylogenetic subtyping using all 491 sequences revealed that sequences fell into one of the three subtype C clades (S1 Fig). Among the 301 sequences from Gondar, 59% (n = 177) were classified as C-EA, 32% (n = 97) as C’-ET and 9% (n = 27) as C-SA. Notably, sequences of the C-SA clade have not been described in Ethiopia previously and the proportion of this clade in Gondar increased from 2% in 2003 to 23% in 2012 (p<0.001, two-tailed linear-by-linear test for association) while there was a modest decline of the two major clades (C-EA: 60%-51%; C’-ET: 38%-26%; Fig 1) during the same time period. However, the proportion of C-EA remained higher than C’-ET and C-SA at each time point 2003–2013. Further analysis of the temporal changes revealed that there was no changes of the proportions of the C-EA and C’-ET clades in Gondar during the same time period (p = 0.146 and p = 0.173, respectively, two-tailed linear-by-linear test for association). Among the sequences with DRMs identified in Gondar (Table 1 and S2 Table), 13 were classified as belonging to the C-EA clade, two to the each of the C’-ET and C-SA clades, however, there were no association of DRMs to either of the clades (C-EA vs C’-ET, p = 0.095; C-EA vs C-SA, p = 1.000; and C’-ET vs C-SA, p = 0.210; FET for all comparisons).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Trend analysis of subtype C clades 2003–2013.

Proportion of different subtype C clades in Gondar 2003–2013. A dotted line indicates the trend over time. The 2012 time point included specimens obtained 2011–2013 (current study).

https://doi.org/10.1371/journal.pone.0205446.g001

Transmission cluster analysis and association with drug resistance mutations

Transmission cluster analysis was performed separately for the C-EA and C-SA/C’-ET sequences by constructing new ML trees of the two data sets (Fig 2). The C-EA phylogeny indicated that the sequences from Gondar as well as sequences from other cities in Ethiopia intermixed with the global dataset (mostly with sequences from other East African countries and Europe). The C-SA/C’-ET phylogeny showed a clear separation of the two clades similar to that seen for the subtyping phylogeny of the whole data set (S1 Fig). While the C-SA clade mostly contained sequences of southern African countries, sequences of Gondar were intermixed in the phylogeny. The C’-ET clade consisted mostly of Ethiopian, followed by sequences obtained in Europe and North America. Gondar transmission clusters were defined as statistically supported phylogenetic clusters dominated by Gondar sequences (see Methods, Table 2 and S3 Table). In total, we identified 28 clusters (which together contained 35% [105 of 301] of the sequences from Gondar) which was indicative for multiple introductions of HIV into Gondar followed by local spread (21 dyads, six medium sized clusters/networks and one large clusters; Fig 2, Table 2 and S3 Table). The Gondar C-EA sequences were found more frequently in transmission clusters compared to the C’-ET sequences (38% [68 of 177 sequences] in 15 local transmission clusters versus 23% [22 out of 97 sequences] in 11 clusters); p = 0.010, FET). Likewise, C-SA sequences were found more frequently in clusters (56% [15 of 27]) compared to C’-ET (p = 0.002, FET) while there were no difference between the number of C-EA and C-SA sequences in cluster (p = 0.098, FET). The mean cluster size for C-EA was five Gondar sequences/cluster (range: 2–28). Nine C-EA clusters had two sequences/cluster (dyads), five had ≥3–8 sequences/cluster while one cluster was large and contained 32 sequences (Fig 2). All C’-ET clusters were represented as dyads while the C-SA clusters was represented by one dyad and one cluster with 13 members. Among C-EA, six of 15 clusters contained >2 sequences which was higher than for C’-ET (none of the 11 clusters had more than two sequences; p = 0.020, FET) but this difference was not apparent between C-EA and C-SA (one cluster contained >2 sequences; p = 1.000, FET). Nine of the 15 C-EA clusters, seven of the 12 C’-ET clusters and one of the two C-SA clusters contained sequences of participants obtained in more than one of the surveys conducted in Gondar 2003–2013 indicated that clusters were populated over time in Gondar (Table 2). Moreover, eight of the 28 clusters (three C-EA clusters, three C’-ET and both of the C-SA clusters) had members from the last survey conducted 2011–2013, indicating recent transmissions within these clusters.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Maximum likelihood phylogenetic tree.

Monophyletic clusters are shown in greyed boxes, defined as clusters with a branch support (aLRT-SH) >0.9 and with >80% sequences obtained from Gondar. Tips are coloured according to collection place. The colour code is indicated in the top of the left panel. A red square at a tip indicates a DRM sequence. (Left panel) The C-SA and C’-ET clades. The branch separating the two clades are indicated with an asterisk and is represented by a branch support of >0.9. (Right panel) The C-EA clade phylogeny.

https://doi.org/10.1371/journal.pone.0205446.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Characteristics of Gondar transmission clusters.

https://doi.org/10.1371/journal.pone.0205446.t002

Dating of transmission clusters of different subtype C clades in Gondar

To study the timing of introductions and differences between the different subtype C clades clusters, we performed Bayesian coalescent analysis using BEAST with sequences obtained from Gondar only (n = 301). Since the Maximum likelihood phylogeny showed that the strains of the different clades in Gondar were dispersed among sequences of other origins and most likely represented multiple introductions, the estimated root tMRCAs of each subtype C clade represent the date of the origin of the circulating subtype C clades in the region (although only Gondar sequences were used in the date estimations) while the estimated tMRCA of the transmission clusters should approximate introductions and local spread of the viral strains in Gondar, cf references [54, 62]. The age of the transmission clusters ranged from estimated dates of 1980–2009. The oldest clusters were estimated to 1980, 1983 and 1990 for C-EA, C-SA and C’-ET, respectively (Table 2 and S1 Table). Fifteen of the 28 clusters had estimated tMRCAs before year 2000. The 15 C-EA and two C-SA transmission clusters had the oldest estimated median tMRCAs (median dates: 1986 and 1991, respectively), while the 11 C’-ET clusters were younger (median date: 2001). Two C’-ET clusters were young (≤5 years) whereas the large C-EA and network C-SA clusters were estimated to be 28 and 30 years with estimated median tMRCAs of 1985 and 1983, respectively. Taken together, this indicated that the C-EA and C-SA was introduced earlier into Gondar compared to the C’-ET strain. The median population evolutionary rate of C-EA was higher than for C’-ET and C-SA (2.29 x10^-3 and 1.85x10^-3 versus 1.70x10^-3 substitutions/site/year, respectively), however the 95% credible intervals overlapped (S1 Table).

Transmission cluster associated with drug resistance mutations

Sequences with DRMs were found more frequently in clusters among the C-EA compared to the C’-ET sequences (57% versus 25%; p<0.001, FET). The C-EA Cluster 12 represented the largest cluster and contained several sub-clusters indicating a large transmission network in Gondar over an extended time period. In subcluster 12a, four of nine sequences carried the G190A mutation (Table 2 and Fig 2) which suggested cluster-associated transmission of this resistance mutation in Gondar. To verify the age estimate of cluster 12 we obtained dated phylogenies using the Bayesian MCMC method by using sequences of this cluster only in the analysis (S3 Table). The median estimated evolutionary rate of the cluster was 1.34×10⁻³ substitutions/site/year, slower than the median rate estimate of the entire C-EA cluster but with overlapping credible intervals (S1 Table). Moreover, the median tMRCA of cluster 12 was estimated to 1980 (95% HPD: 1969–1989) which was similar to the estimate (1985; 95% HPD: 1967–1985) obtained from the analysis of all Gondar sequences (Table 2 and S1 Table). Since the internode intervals in the MCMC tree represent the maximum transmission time intervals, we estimated the time between transmissions events. Analysis of the internode intervals (from the MCMC tree; S2 Fig) revealed that the maximum median transmission intervals for cluster 12 were 1.58 years (IQR: 1.24–2.31). The median parametric estimate of the population growth rate also indicated frequent transmissions associated with sequences of cluster 12 at approximately 0.5 new infections/individual/year which correspond to a median doubling time of infections every 1.34 years (16 months). The demographic plot analysis indicated that this cluster increased in size and the number of effective infections, N_e, reached its peak in 2003, but declined thereafter (from 4990 in 2003 to 2404 in 2013; S1 Table). However, the four sequences that had the G190A mutation were collected during three of the four surveillance studies in Gondar 2003–2013 (Table 2 and Fig 2) suggesting onward transmission of this DRM over many years in Gondar.