New sequences were generated from blood plasma obtained through studies conducted through the MSM Health Research Consortium—a multisite collaboration between researchers affiliated with KEMRI-Wellcome Trust (KWTRP) in Coastal Kenya, Nyanza Reproductive Health Society (NRHS) in Nyanza, and Sex Workers Outreach Program (SWOP) clinics in Nairobi. These included samples from Coast derived from participants in a prospective observational cohort (2006–2019) (Sanders et al., 2013), samples from Nairobi from a respondent-driven sample survey (Transform, 2017; Smith et al., 2021b), and samples from Nyanza derived from the Anza Mapema cohort (2015–2017) (Kunzweiler et al., 2018).

The HIV-1 pol sequences were comprised of 1,020 nucleotides, HXB2 [K03455] positions 2267–3287. HIV-1 RNA was purified from patient blood plasma using the RNeasy Lipid Tissue Mini Kit (QIAGEN) as previously described (Esbjörnsson et al., 2010). Reverse transcription and amplification of partial pol gene were performed using the One-Step Superscript III RT/Platinum Taq High Fidelity Enzyme Mix (Thermo Fisher Scientific™) with the pol-specific primer pair JA269 and JA272 (Hedskog et al., 2010). First-round PCR products were amplified in a nested PCR with DreamTaq Green DNA Polymerase (Thermo Fisher Scientific™) using pol-specific primers JA271 and JA270 (Hedskog et al., 2010). PCR products were sequenced in both directions with the nested PCR primers using the BigDye terminator kit v1.1 (Applied Biosystems), and the sequences were determined on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems).

Additional Kenyan HIV-1 pol sequences (referred to as published sequences, 2006–2019) were retrieved (October 11th 2021) from the Los Alamos HIV-1 sequence database (Los Alamos National Laboratory, 2019). The combined new and published sequences (referred to as the Kenyan dataset) were annotated with information on sampling dates and geographical area of residence during sampling (i.e., province; Coast, Nairobi, and Nyanza).

The Kenyan dataset was aligned with the HIV-1 Group M (subtypes A-K + Recombinants) subtype reference dataset (available at the Los Alamos HIV database)¹ using the MAFFT algorithm in Geneious Prime 2019 (Larkin et al., 2007; Los Alamos National Laboratory, 2019). The resulting alignment was used to construct a maximum-likelihood phylogenetic tree in PhyML using the general time-reversible substitution model with a gamma-distributed rate variation and proportion of invariant sites (GTR+Γ4+I) (Guindon et al., 2005). Branch support was assessed using the Shimodaira–Hasegawa like approximate likelihood ratio test (aLRT-SH) in PhyML, with aLRT-SH ≥ 0.90 considered as significant (Anisimova and Gascuel, 2006; Hassan et al., 2017). Subtypes were assigned based on the Subtype/CRF-resolved phylogeny visualized using FigTree v1.4.4². Subtype assignment was further verified using the REGA HIV-1 Subtyping Tool (v.3.0), and unique recombinant forms (URFs) were detected using the jumping profile hidden Markov model (jpHMM) (Schultz et al., 2009; Pineda-Peña et al., 2013).

Sequences were grouped into subtype-specific datasets, and a search for related sequences was done for each subtype-specific (A1, C, and D) dataset using the NCBI GenBank BLAST tool, limiting results to the 10 most similar hits per sequence, and retaining the oldest sequence per individual (Kouyos et al., 2010; Esbjörnsson et al., 2016; Sallam et al., 2017). Kenyan sequences and reference sequences were combined and aligned using the MAFFT algorithm in Geneious Prime 2019 (Larkin et al., 2007). Subtype-specific alignments were edited to exclude codon positions associated with drug resistance, and maximum-likelihood phylogenies were reconstructed in PhyML. For each subtype, monophyletic clades with aLRT-SH support ≥0.9 and which were dominated (≥80%) by Kenyan sequences (compared to reference sequences) were defined as Kenyan HIV-1 clusters (Hassan et al., 2017). Clusters were classified based on the number of sequences per cluster into dyads (2 sequences), networks (3–14 sequences), and large clusters (>14 sequences) (Esbjörnsson et al., 2016; Nduva et al., 2020; Abidi et al., 2021).

To date clusters and to estimate the effective population size through time (N_e.T), Bayesian phylodynamic inference was performed in BEAST 1.10.4 using the Bayesian Skygrid model, an uncorrelated lognormal relaxed clock, and the general time-reversible substitution model with a gamma-distributed rate variation and proportion of invariant sites (GTR+Γ4+I) (Drummond et al., 2005; Baele et al., 2012; Gill et al., 2013; Suchard et al., 2018). Only sequences classified as pure A1, C, and D subtypes were analyzed. BEAST runs were computed with a chain length of 100–300 million generations for each dataset, sampling every 10,000th–30,000th iteration, and discarding the first 10% as burn-in. Convergence was determined in Tracer v.1.7.0 and defined as effective sample sizes (ESS) ≥100 (Suchard et al., 2018). Maximum clade credibility (MCC) trees were summarized using Tree-Annotator v1.8.2 (BEAST suite).

To infer the direction of virus movements between geographic locations from HIV-1 sequence data, a discrete phylogeographic inference was computed using specific locations as independent discrete states (Lemey et al., 2009; Edwards et al., 2011; Faria et al., 2014). Several sensitivity analyses were performed to test the robustness of our data. Firstly, the Kenyan dataset was grouped by subtype (A1, C, and D), and the phylogeographic inference was performed using all the sequences per subtype. Secondly, to reduce sampling bias arising from the unproportionable allocation of sequences per location, sequences in the subtype A1-specific dataset (the largest of the three subtypes) were randomized and subsampled into a dataset with an equal number of sequences per province using in-house Perl scripts (available upon request). Lastly, subtype A1 sequences from Coast province were subsampled uniformly and used to estimate virus migration between three geographically distinct regions in Coastal Kenya (i.e., Mombasa, South Coast, and North Coast).

In the phylogeographic inference, the asymmetric model was adopted (over the alternative symmetric model) as it relaxes the assumption of constant diffusion rates through time to realistically model the location-exchange processes (Edwards et al., 2011; Faria et al., 2014). In addition to estimating the direction of HIV-1 migration, the proportions of forward and reverse rates of migrations between geographic locations were quantified using a robust counting approach (Markov jumps) implemented in BEAST (Minin and Suchard, 2008). Maximum clade credibility (MCC) trees annotated with demographic and epidemiological data were summarized in Tree-Annotator v1.10.4 (BEAST suite) and visualized in FigTree (v1.4.4). Well-supported virus movements and Bayes factors (BFs) assessing statistical support were summarized using SPREAD v1.0.7, and BF ≥ 3 was considered significant (Lemey et al., 2009).

Continuous data were presented using medians and interquartile ranges (IQRs). Frequencies and percentages were used to describe categorical data. A multivariable logistic regression model was used to assess associations between individual sequence characteristics [e.g., subtype, location of sampling, year (range) of sampling, and source of sequence data—i.e., published or newly generated] and phylogenetic clustering. Statistics and summary plots were done using Stata 15 (StataCorp LLC, College Station, TX, United States) and RStudio (version 1.2.5001) with the packages yarrr and ggplot2 (Wickham, 2016; Phillips, 2017).

Plasma samples used to generate the new sequences were obtained from ongoing or concluded studies that were also approved by Kenya Medical Research Institute (KEMRI) Scientific and Ethics Review Unit (SERU 3747, 3280, and 3520, and SSC 894). Since published sequences were obtained from an open-access public domain, informed consent was not retrospectively obtained. Instead, we sought approval through a study protocol that was reviewed by KEMRI/SERU (SERU 3547).

Among the 372 HIV-1 partial pol sequences analyzed, 213 (57.3%) were generated in this study, and 159 (42.7%) were previously published. The majority (N = 178, 47.9%) of the sequences were from Coast province, 137 (36.8%) from Nairobi province, and 57 (15.3%) from Nyanza province (Figure 1, Table 1, Supplementary Figure 1, and Supplementary Tables 1, 2). Sequences belonged to sub-subtype A1 (N = 268, 72.0%), subtype D (N = 41, 11.0%), subtype C (N = 22, 5.9%), subtype G (N = 3, 0.8%), CRF 21A2D (N = 3, 0.8%), CRF 16A2D (N = 1, 0.3%), and subtype B (N = 1, 0.3%). Unique recombinant forms (URFs) identified included A1D (N = 19, 5.1%), A1C (N = 7, 1.9%), D01AE (N = 5, 1.3%), A1B (N = 1, 0.3%), and DB (N = 1, 0.3%, Figure 2).

Clusters were determined from maximum-likelihood (ML) phylogenies reconstructed for the most prevalent HIV-1 subtypes in the population [subtypes A (A1), C, and D—cumulatively comprising 89.0% of the sequences in the Kenyan dataset]. Non-Kenyan HIV-1 reference sequences were obtained from GenBank based on similarity (where of 931 participant-unique sub-subtype A1 sequences remained after removal of redundancies; 488 for subtype C; and 350 for subtype D). Of 331 (A1, C, and D) sequences in the cluster analysis, 229 sequences (61.2%) formed 46 statistically supported clusters (size range: 2–20 sequences). Dyad/pairs were most common (N = 25, 54.4% of all clusters), followed by networks having 3–14 sequences (N = 18, 39.1%), and large clusters having more than 14 sequences (N = 3, 6.5%). The majority (N = 34, 73.9%) were sub-subtype A1 clusters, followed by subtype D (N = 8, 17.4%) and subtype C (N = 4, 8.7%, Table 2 and Supplementary Figure 2).

Stratification of clusters by geographic regions showed two distinct clustering patterns. First, some clusters (N = 23, 50.0%) had sequences belonging exclusively to one specific province including Coast (N = 14, 30.4%), Nairobi (N = 6, 13.0%), and Nyanza (N = 3, 6.5%) province-exclusive clusters. The remaining clusters (N = 23, 50.0%) were mixed between different provinces where HIV-1 mixing between Coast and Nairobi was most common (N = 13, 28.3% clusters), followed by mixing between Nyanza, Nairobi, and Coast (N = 5, 10.9%), Nyanza and Nairobi (N = 3, 6.5%), and Nyanza and Coast (N = 2, 4.4%, Table 2 and Supplementary Figure 2). Sequences from Nairobi province were more likely to cluster compared to sequences from Coast province [adjusted odds ratio (aOR) 3.5, 95% confidence interval (CI) 1.2–10.4, P = 0.022, Table 3].

In-depth phylodynamic analysis indicated that the number of MSM contributing to new HIV-1 A1 infections over time increased exponentially during the early 2000s, followed by a period with some fluctuation (but largely steady) between 2000 and 2017, and mostly decreasing dynamics during recent years (2017–2019, Figure 3A). Likewise, for both subtype C and D lineages, the effective population size increased exponentially during 2007–2008 and has stabilized in recent years (2016–2019; Figures 3B–D).

Estimating dates of origins of all clusters indicated that the majority (65%) of transmissions within clusters took place between 2000 and 2014. The oldest sub-subtype A1 cluster had nine MSM from Nyanza, Nairobi, and Coast and had originated during 1987, while the youngest cluster was dated to 2014 among MSM in Nyanza (Figure 4A, Supplementary Figure 3, and Supplementary Table 3). The largest A1 cluster (N = 20, 2008–2017) had remained active over 20 years since the estimated time to the most recent common ancestor (tMRCA) in 1997 and was geographically spread out to Nyanza, Nairobi, and Coast provinces. The second-largest A1 cluster (N = 19, 2008–2017) originated in 1996 and had sequences from Nyanza, Nairobi, and Coast provinces. The four subtype C clusters originated during 1988, 1998, 2009, and 2014, respectively, while the earliest subtype D cluster originated during 1976 and the youngest during 2014 (Figures 4B,C and Supplementary Table 3). Overall, there was evidence of onward HIV-1 transmission among MSM, within longstanding and geographically diverse HIV-1 networks.

Ancestral locations and rates in historical virus jumps were first estimated based on all subtype-specific sequences in the Kenyan dataset (i.e., 268 sub-subtype A1, 41 subtype D, and 22 subtype C sequences). Phylogeographic analysis indicated significant support (Bayes Factor, BF ≥ 3) for virus migration from Coast to Nairobi (BF = 3716; subtype A1, BF = 268; subtype C; and BF = 16; subtype D) and from Nairobi to Nyanza (BF = 3716; subtype A1, BF = 43; subtype D, Supplementary Table 4). Exploring temporal trends in virus transitions between geographic provinces summarized from trait-annotated maximum clade credibility trees indicated that the proportion of virus export from Coast to Nairobi increased from 4.2% before 2000 to 14.2% during 2001–2010 and declined to 4.9% during 2011–2020. Likewise, virus export from Nairobi to Nyanza increased from 2.4% in 2000–2010 to 10.8% in 2011–2020, while reverse transitions were rare and occurred only from Nyanza to Nairobi (Supplementary Figures 4, 5 and Supplementary Table 5).

A sensitivity analysis with uniform sampling per province was performed to confirm the robustness of the initial phylogeographic inference. The uniformly subsampled dataset comprised 135 HIV-1 sub-subtype A1 sequences (45 sequences each from Nairobi, Mombasa, and Nyanza province). Based on this analysis, there was significant support for HIV-1 migration from Coast to Nairobi (BF = 7766), Nairobi to Nyanza (BF = 1293), and Coast to Nyanza (BF = 336, Table 4). Furthermore, Markov jumps estimates with uniform sampling indicated that the majority (80.3%) of HIV-1 jumps between provinces occurred from Coast to other provinces including jumps from Coast to Nyanza (N = 26, 42.6% of all virus jumps between provinces) and from Coast to Nairobi (N = 23, 37.7%, Figure 5 and Table 5). There was also some (N = 10, 16.4%) virus exchange between Nairobi and Nyanza, such that virus jump Nairobi to Nyanza (N = 7, 11.5%) was twofold higher than from Nyanza to Nairobi (N = 3, 4.9%, Table 5).