Temporal trends and transmission dynamics of pre-treatment HIV-1 drug resistance within and between risk groups in Kenya, 1986–2020

Abstract Background Evidence on the distribution of pre-treatment HIV-1 drug resistance (HIVDR) among risk groups is limited in Africa. We assessed the prevalence, trends and transmission dynamics of pre-treatment HIVDR within and between MSM, people who inject drugs (PWID), female sex workers (FSWs), heterosexuals (HETs) and perinatally infected children in Kenya. Methods HIV-1 partial pol sequences from antiretroviral-naive individuals collected from multiple sources between 1986 and 2020 were used. Pre-treatment reverse transcriptase inhibitor (RTI), PI and integrase inhibitor (INSTI) mutations were assessed using the Stanford HIVDR database. Phylogenetic methods were used to determine and date transmission clusters. Results Of 3567 sequences analysed, 550 (15.4%, 95% CI: 14.2–16.6) had at least one pre-treatment HIVDR mutation, which was most prevalent amongst children (41.3%), followed by PWID (31.0%), MSM (19.9%), FSWs (15.1%) and HETs (13.9%). Overall, pre-treatment HIVDR increased consistently, from 6.9% (before 2005) to 24.2% (2016–20). Among HETs, pre-treatment HIVDR increased from 6.6% (before 2005) to 20.2% (2011–15), but dropped to 6.5% (2016–20). Additionally, 32 clusters with shared pre-treatment HIVDR mutations were identified. The majority of clusters had R0 ≥ 1.0, indicating ongoing transmissions. The largest was a K103N cluster involving 16 MSM sequences sampled between 2010 and 2017, with an estimated time to the most recent common ancestor (tMRCA) of 2005 [95% higher posterior density (HPD), 2000–08], indicating propagation over 12 years. Conclusions Compared to HETs, children and key populations had higher levels of pre-treatment HIVDR. Introduction of INSTIs after 2017 may have abrogated the increase in pre-treatment RTI mutations, albeit in the HET population only. Taken together, our findings underscore the need for targeted efforts towards equitable access to ART for children and key populations in Kenya.


Introduction
By the end of 2020, an estimated 30 million individuals were receiving ART globally. 12][3] However, widespread use of ART has been associated with the emergence of HIV-1 drug resistance (HIVDR) and the transmission of drug-resistant viruses that can compromise therapy outcomes. 4he WHO recommends routine nationally representative HIVDR surveys to inform treatment strategies. 5Increasing levels of pretreatment HIVDR to NNRTIs has been recognized as a problem for many years, especially in low-and middle-income countries (LMICs) where virological monitoring of patients on ART is suboptimal. 6,7In part, the increasing NNRTI resistance and poor tolerability to PIs informed WHO's recommendation for a switch to integrase strand transfer inhibitors (INSTIs) as part of the first-line regimens in LMICs.INSTIs are more efficacious and have a high genetic barrier to resistance.Dolutegravir, a second-generation INSTI, is now widely adopted for both treatment-naive and experienced PWHIV. 8owever, the potential impact of the switch to dolutegravir-based regimens, which came into effect in LMICs in 2017, on emergence and circulation of NNRTI-resistant strains is unclear.
Further, key populations including MSM, people who inject drugs (PWID), female sex workers (FSWs) and transgender people are disproportionately impacted by HIV-1 in sub-Saharan Africa (sSA).HIV-1 transmission from these key populations is thought to contribute disproportionately to Africa's generalized HIV-1 epidemic. 9e previously demonstrated that the majority of HIV-1 transmissions in Kenya occur within distinct risk groups, and that transmission from heterosexuals (HETs) to key populations is more common than vice versa. 10,11However, it remains unknown whether there is differential transmission of HIVDR mutations among different risk groups.We aimed to describe prevalence, temporal trends and transmission linkages of pre-treatment HIVDR mutations within and between risk groups in Kenya.

Study population
HIV-1 partial pol sequences from Kenya were either newly generated from archived plasma samples (henceforth referred as newly generated) or retrieved from the Los Alamos HIV-1 sequence database (henceforth referred as published sequences). 12Newly generated HIV-1 pol sequences [approximately 1020 nucleotides (nt), HXB2 (K03455) positions 2267-3287] were generated from plasma obtained from the MSM Health Research Consortium.This is a multi-site collaboration between the Kenya Medical Research Institute-Wellcome Trust Research Programme (KWTRP) in Coastal Kenya, [13][14][15] Nyanza Reproductive Health Society (NRHS) in Western Kenya, 16 Sex Workers Outreach Program (SWOP) clinics in Nairobi, Kenya, the Targeted Research Advancing Sexual Health for Men who have Sex with Men (TRANSFORM) cohort from Nairobi, Kenya, 17 and the national HIV-1 reference laboratory at the Kenya Medical Research Institute/Centre for Global Research (KEMRI/CGHR) in Kisumu, Kenya.Plasma samples from the collaboration were used to generate HIV-1 partial pol sequences as previously described. 18Published HIV-1 pol sequences of Kenyan origin were also retrieved from the Los Alamos HIV database on 20 March 2022. 12ociodemographic, clinical and ART status data for newly generated sequences were obtained from the respective cohorts, while those from published sequences were extracted from referenced publications and their accompanying supplementary material.Where data for published sequences were missing, information was obtained through direct communication with study authors.Newly generated and published sequences were then annotated with sampling date, sampling location, risk group [HETs, MSM, PWID, FSWs and children (defined as individuals <18 years old)] and ART status (whether treatment naive or experienced).Where such information was still missing, related sequences were excluded from the analysis.Only sequences generated from ART-naive individuals were included in the analyses.

HIV-1 subtype determination
Newly generated and published HIV-1 partial pol sequences were aligned with the HIV-1 Group M (subtypes A-K + recombinants) subtype reference sequences (http://www.hiv.lanl.gov)using the MAFFT algorithm in Geneious Prime 2019. 12Subtypes were determined based on a maximum-likelihood (ML) phylogenetic tree generated in PhyML using the general time-reversible substitution model with a gamma-distributed rate variation and proportion of invariant sites (GTR + Γ4 + Ι). 19 Branch support was determined using the approximate likelihood ratio test with the Shimodaira-Hasegawa-like procedure (aLRT-SH) in PhyML and aLRT-SH ≥ 0.90 was considered significant.The phylogenetic tree was visualized using FigTree v1.4.4 (https://github.com/rambaut/figtree/releases),and subtypes assigned based on how sequences in our dataset clustered relative to the HIV-1 Group M reference sequences.Subtype assignment congruence was confirmed with the web-based subtyping algorithm COMET (available at https://comet.lih.lu) and REGA (version 3.46, https://www.genomedetective.com/app/typingtool/hiv).[22]

Characterization of pre-treatment HIVDR
All sequences from antiretroviral-naive individuals were submitted to the Stanford HIV database (https://hivdb.stanford.edu)for determination and interpretation of HIVDR mutations.The Calibrated Population Resistance tool for protease and reverse transcriptase (PRRT) and for integrase (IN) sequences was used to identify pre-treatment HIVDR mutations based on the WHO list of mutations for the surveillance of transmitted HIVDR. 23,24Mutations were grouped and presented based on drug class as follows: NRTI, NNRTI, PI and INSTI mutations.

Phylogenetic determination of clusters
HIV-1 sequences were grouped into subtype-specific datasets.For each index sequence, the 20 most similar sequences were obtained from the NCBI GenBank using the BLAST tool, as previously described. 25After removal of duplicate sequences from the same individual, the Kenyan sequences from our dataset and the reference sequences obtained from the BLAST search were compiled into subtype-specific alignments, excluding codon positions associated with drug-resistance mutations.ML phylogenies were generated in PhyML using the GTR + Γ4 + Ι model. 19For each subtype-specific phylogenetic tree, monophyletic clades with aLRT-SH support ≥ 0.9 and which were dominated (≥80%) by Kenyan sequences (compared with reference sequences) were defined as Kenyan HIV-1 clusters.Clusters were classified based on the number of sequences into dyads (two sequences), networks (3-14) and large clusters (>14 sequences). 25ML phylogenetic trees were annotated based on the most common mutations per drug class, i.e.K103NS, Y181CIV and G190AES (for NNRTI resistance mutations) and M184VI and T215 revs (for NRTI resistance mutations).

Hierarchical phylogenetic modelling and Bayesian inference
Hierarchical phylogenetic models (HPMs) are useful for phylogenetic inferences when analysing a sparse sample. 26HPMs comprise two components: (i) an across-partition-level model, which captures the shared evolutionary history among all species in an analysis, often at a higher taxonomic level; and (ii) an individual partition model that provides more detailed information about evolutionary relationships within specific partitions.HPMs innovatively allow for the simultaneous fitting of both the across-partition-level model and individual partition models.Thus, information inferred at the broader taxonomic level can be shared as a prior with the individual partition models.This 'feedback' or sharing of information results in a borrowing Nduva et al.  of strength from one partition to another, which can lead to more accurate and precise estimates of evolutionary relationships. 26PMs were incorporated in our phylodynamic analyses to date clusters.BEAST software (version 1.10.4) was used.The HKY nucleotide substitution model (unlinked for all partitions) was applied while selecting empirical base frequencies and a gamma site heterogeneity model with four gamma categories.Codon partitions were turned off.A strict clock model (unlinked for all partitions) was specified, and a constant size tree prior (with unlinked tree priors for all trees) maintained as a simple coalescent tree prior.We specified hierarchical priors for the substitution model parameters (on the log-normal scale for both 'kappa' and 'alpha' parameters, while setting the normal population mean hyperprior standard deviation to 5.0, as well as a gamma population precision hyperprior of initial value 1.0, shape value 0.01 and scale value 100.0).A second HPM was set for the clock rates (on the log-normal scale, with a normal population mean hyperprior standard deviation set to 100.0, and a gamma population precision hyperprior with an initial value 1.0, shape value 0.001 and scale value 1000.0).A final HPM was specified for the population sizes, setting the normal hyperprior standard deviation to 100.0, and the gamma hyperprior shape and scale to 0.001 and 1000.0,respectively.Markov chain Monte Carlo (MCMC) runs with a chain frequency of 300 million generations were computed in BEAST, logging every 300 000th iteration, and discarding the first 10% as chain burn-in.Convergence was determined in Tracer v.1.7.0 [defined as effective sample size (ESS) ≥ 100]. 27Maximum clade credibility (MCC) trees were summarized from the posterior tree distribution using TreeAnnotator v1.10.4 (BEAST suite) and visualized using Figtree (v1.4.4,https://github.com/rambaut/figtree/releases).
9][30] In the absence of empirical data on the duration of infectiousness in Kenya, R 0 was modelled assuming varying D values (range, 1-8 years).

Statistical analysis
Prevalence of pre-treatment HIVDR was determined as a proportion of sequences with at least one HIVDR mutation relative to the total number of sequences included in the analysis.Prevalence estimates were aggregated by calendar year of sampling into 5 year intervals to reflect major changes in treatment guidelines in Kenya as follows: before 2005 (before introduction of combination ART), 2005-10 (introduction of combination ART), 2011-15 (scale up of combination ART) and 2016-20 (introduction of INSTI-based regimens).Categorical data were compared with the χ² test and continuous data were compared with the Kruskal-Wallis test, where appropriate.Statistical differences between risk group and calendar year intervals were assessed using the Dunn's test for multiple comparisons, with Bonferroni correction.Overall, drug class-specific, and mutation-specific HIVDR prevalence estimates and 95% CIs were presented.HIVDR temporal trends were assessed using nptrends, a non-parametric extension of the Wilcoxon rank-sum test. 31Data analysis and summary plots were done using Stata 15 (StataCorp LLC, College Station, TX, USA) and RStudio (version 1.2.5001) with the ggplot2 package. 32

Nucleotide sequence accession numbers
Newly generated nucleotide sequences are available from GenBank under the accession numbers (MT084914-MT085076) and (OM109695-OM110282).

Ethics
For newly generated sequences, informed consent for use of data and samples for research studies was obtained from participants under respective study protocols.Since published sequences were obtained from an open access repository, informed consent was not retrospectively obtained.Instead, science and ethics approval were obtained from the Kenya Medical Research Institute (KEMRI) Scientific and Ethics Review Unit (SERU 3547).All newly generated and published sequences were de-identified at source before inclusion in the analysis.
Bayesian dating was performed for networks and large clusters, which included four sub-subtype A1 and one subtype D clusters (Table 3 and Figure 5).The largest comprised 16 MSM with the K103N mutation.This cluster had an estimated time to the most recent common ancestor (tMRCA) of 2005 [95% higher posterior density (HPD), 2000-08], where the most recent sample was collected in 2017, suggesting that this lineage had persisted   3 and Table S8).

Discussion
We used HIV-1 pol sequence data spanning a period of more than three decades to assess the prevalence, temporal trends and transmission of pre-treatment HIVDR among different risk groups in Kenya.Overall prevalence of pre-treatment HIVDR was high (15.4%),which was largely a reflection of the high levels of pre-treatment NNRTI drug-resistance mutations.Notably, there was comparatively lower (<10%) pretreatment NRTI drug resistance and no pre-treatment INSTI resistance in the study population, albeit based on a limited number of HIV-1 IN gene sequences mostly collected prior to the transition to INSTI-based regimens in Kenya.Taken together, these observations justify the retention of two NRTIs and the switch to a dolutegravir-based regimen for HIV-1 treatment in Kenya. 33e observed an increasing trend in pre-treatment HIVDR at a countrywide level, which was also largely driven by NNRTI resistance.][36] NNRTIs have a low genetic barrier to resistance, which permits outgrowth of resistant variants under suboptimal drug pressure, and mutations may persist for long durations, facilitating their onward transmission. 37,38Of interest, pre-treatment NNRTI and NRTI resistance among HETs (the largest risk group in this study, which is representative of the HIV-1 epidemic in Kenya) increased between 2005 and 2015 but declined between 2015 and 2020.This coincides with the nationwide transition from NNRTI to INSTI-based ART regimens, 33 suggesting that scale-up of dolutegravir-based regimens may have abrogated emergence and circulation of pre-treatment NNRTI and NRTI resistance mutations.In contrast, pre-treatment HIVDR among FSWs and MSM increased consistently through to 2015-20.Interestingly, as of 2020, ART coverage was lower among key populations: 73% in  HET is defined as presumed heterosexual, i.e. at-risk men and women not reporting sex work or male same-sex behaviour.Pre-treatment HIV-1 drug resistance in Kenya (1986-2020)  Nduva et al.
FSWs, 68% in PWID and 63% in MSM, compared with 86% in the general HET population. 39Taken together, these observations may explain the higher resistance levels observed in key populations compared with the HET population, and warrants further investigation.
Our observations are consistent with findings from a WHO-led systematic review reporting higher levels of pre-treatment HIVDR among key populations compared with the HET population in the global context (with better sampling densities from higherincome countries, but with limited representation of data from key populations in Africa). 5With the transition to INSTI-based regimens, enhanced monitoring of INSTI resistance using contemporary sequences from children, HETs and key populations would be prudent to inform future treatment strategies in Kenya.
We observed high levels of K103NS, Y181CIV and G190AES mutations, which likely reflect extensive selection in persons receiving nevirapine and efavirenz, both of which comprised the main NNRTI options for first-line regimens historically in Kenya. 40This observation was more evident among children, who had the highest prevalence of pre-treatment HIVDR, with NNRTI-associated mutations Y181CIV and K103NS being the most common, which is consistent with literature. 41,42These mutations likely point towards the widespread use of nevirapine for prevention of mother-to-child transmission (pMTCT).Due to a paucity of HIV-1 sequence data from children during the 2016-20 time interval, temporal trend analyses were not possible.Continued surveillance of pre-treatment HIVDR in children to assess the impact, if any, of the transition to a dolutegravir-based regimen on NNRTI-associated mutations in Kenya is therefore warranted.
We observed several HIV-1 clusters with shared mutations spanning more than 10 years, indicating onward propagation of HIVDR among treatment-naive individuals as new infections, which has also been reported from developed settings. 43There was more frequent HIV-1 clustering of K103NS strains, indicating extensive transmission of this pre-treatment mutation in Kenya, possibly dating back to the single-dose nevirapine era.Most of the clusters had a basic reproductive number, R 0 ≥ 1.0, suggesting ongoing propagation of NRTI and NNRTI HIVDR transmissions.Assessment of these trends over the next years is needed to determine whether roll-out of the more efficacious INSTIs as first-line regimens will terminate circulation of NNRTI resistance mutations in Kenya.
The main strength of our study was the analysis of sequence data from multiple geographically diverse regions from Kenya and spanning over three decades.However, the study is not without limitations, most of which are consistent with populationlevel studies that leverage secondary sources of data.First, although we included all the Kenyan HIV-1 pol sequences available in the public domain and complemented the sample size by generating more sequences from well-characterized cohorts, the distribution of sequences was still sparse in some risk groups during some time intervals.This may have resulted in potential sampling bias, which may have impacted our temporal trends and phylodynamic inferences, and warrants cautious interpretation of our findings.Second, we extracted a limited set of demographic and clinical data including year of sampling, risk group and ART status.This implies that we were not able to control for other factors like duration from the estimated date of infection and unreported history of preexposure prophylaxis use.While it would have been ideal to control for these factors in our analyses, most are usually either not available or not systematically collected across parent studies.
In conclusion, we demonstrated an increase in pre-treatment HIVDR over the last two decades, which justifies the switch to INSTI-based therapy.Importantly, we show that children and key populations have significantly higher levels of pre-treatment HVIDR compared with the general HET population in Kenya.Remarkably, INSTI-based regimens may have abrogated propagation of RTI mutations among the general HET population, but Pre-treatment HIV-1 drug resistance in Kenya (1986-2020)  not in hard-to-reach key populations.Our study further underscores the challenges with access to HIV care and treatment for key populations, as demonstrated with increasing pretreatment HIVDR in MSM and FSWs, despite the broad nationwide introduction of INSTI-based regimens.Using phylodynamic analyses, we also demonstrate long-standing and ongoing propagation of pre-treatment HIVDR strains, that were mostly risk-group exclusive.Taken together, our findings underscore need for targeted efforts towards equitable access to ART for children and key populations in Kenya.Bayesian dating was restricted to five clusters having ≥3 sequences and with the three most dominant pre-treatment HIVDR mutations per drug class, i.e.NNRTI (K103N and Y181C) and NRTI mutations (T215revs).One cluster shared more than one mutation (K103N/K101P).HET is defined as at-risk men and women who did not report sex work or male same-sex behaviour.Nduva et al.

Figure 1 .
Figure 1.A flowchart summarizing the number of HIV-1 pol sequences(1986-2020) analysed in this study.RTI, reverse transcriptase inhibitor.The HET population was defined as at-risk men and women who did not report sex work or male same-sex behaviour.

Figure 2 .
Figure 2. HIV-1 subtypes distribution within and between risk groups in Kenya (n = 3567; 1986-2020).A stacked bar graph summarizing subtype distribution in the full Kenyan HIV-1 sequence dataset.Subtypes are depicted on the x-axis and the frequency per risk group on the y-axis, coloured per risk group (green: children; sky blue: FSWs; vermillion: HETs; purple: MSM and dark blue: PWID).The most dominant subtypes were A1, D and C (and their recombinant forms) (n = 3567; 1986-2020).

Figure 3 .
Figure 3. Distribution of pre-treatment HIVDR by: (a) proportion of individuals with at least one pre-treatment HIVDR mutation by risk groups, with pairwise comparisons and Bonferroni corrections; (b) proportion of individuals with at least one pre-treatment HIVDR mutation by drug class (with 95% CI); (c) proportion of individuals with at least one pre-treatment HIVDR mutation by both risk group and drug class (with 95% CI); (d) proportion of pre-treatment HIVDR mutations by drug class (with 95% CI); and (e) proportion of the five most prevalent pre-treatment HIVDR mutations per drug class distributed by risk groups, in Kenya (n = 3567; 1986-2020).

Figure 4 .
Figure 4. Temporal trends in pre-treatment HIVDR in Kenya among (a) overall ART-naive individuals (all risk groups combined); (b) HETs in the general population; (c) FSWs; and (d) MSM.Proportions are presented with error bars representing 95% CIs as estimated by a non-parametric linear-by-linear test for trends (nptrends) (n = 3567; 1986-2020).

Figure 5 .
Figure 5. MCC trees estimated using a Bayesian MCMC approach showing Kenyan HIV clusters (≥3 sequences) with shared pre-treatment HIVDR mutations.For each cluster, the estimated time of origin (95% HPD intervals) is indicated at the root node, and branch tips correspond to sampling dates per sequence.Legends correspond to cluster number (#), subtype (A1 or D), mutation and risk group.

Table 1 .
Distribution of newly generated and published HIV-1 pol sequences from antiretroviral therapy-naive individuals in Kenya (n = 3567; 1986-2020)

Table 2 .
Characteristics of clusters (n = 32) with shared pre-treatment HIVDR mutations and distributed into subtypes and risk group Number of clusters with the most dominant pre-treatment HIVDR mutations per drug class i.e.NNRTI (K103NS, Y181CIV, G190AES) and NRTI (M184VI and T215revs).Two clusters shared more than one mutation (K103NS/M184VI and K103N/K101P).HET is defined as presumed heterosexual, i.e. at-risk men and women not reporting sex work or male samesex behaviour.