Characterization of HIV-1 Epidemic in Kyrgyzstan

Kyrgyzstan has one of the highest rates of HIV-1 spread in Central Asia. In this study, we used molecular–epidemiological approaches to examine the HIV-1 epidemic in Kyrgyzstan. Samples were obtained from HIV-positive individuals who visited HIV/AIDS clinics. Partial pol gene sequences were used to identify HIV-1 subtypes and drug resistance mutations (DRMs) and to perform phylogenetic analysis. Genetic diversity and history reconstruction of the major HIV-1 subtypes were explored using BEAST. This study includes an analysis of 555 HIV-positive individuals. The study population was equally represented by men and women aged 1–72 years. Heterosexual transmission was the most frequent, followed by nosocomial infection. Men were more likely to acquire HIV-1 during injection drug use and while getting clinical services, while women were more likely to be infected through sexual contacts (p < 0.01). Heterosexual transmission was the more prevalent among individuals 25–49 years old; individuals over 49 years old were more likely to be persons who inject drugs (PWID). The major HIV-1 variants were CRF02_AG, CRF63_02A, and sub-subtype A6. Major DRMs were detected in 26.9% of the study individuals; 62.2% of those had DRMs to at least two antiretroviral (ARV) drug classes. Phylogenetic analysis revealed a well-defined structure of CRF02_AG, indicating locally evolving sub-epidemics. The lack of well-defined phylogenetic structure was observed for sub-subtype A6. The estimated origin date of CRF02_AG was January 1997; CRF63_02A, April 2004; and A6, June 1995. A rapid evolutionary dynamic of CRF02_AG and A6 among Kyrgyz population since the mid-1990s was observed. We observed the high levels of HIV-1 genetic diversity and drug resistance in the study population. Complex patterns of HIV-1 phylogenetics in Kyrgyzstan were found. This study highlights the importance of molecular–epidemiological analysis for HIV-1 surveillance and treatment implementation to reduce new HIV-1 infections.

Introduction

Kyrgyzstan, or Kyrgyz Republic, is a small country in Central Asia bordering with Kazakhstan, Uzbekistan, Tajikistan, and China. In 2020, an estimated 9,200 (8,400--9,900) people lived with HIV-1 in Kyrgyzstan with a prevalence of 0.2%¹. The HIV-1 epidemic concentrates in the key populations, mainly persons who inject drugs (PWID) and their sexual partners (Zhao et al., 2020). Access to HIV-1 care and treatment services is far behind the global target of 95-95-95 by 2030. In 2020, 76% of people were aware of their HIV-1 status, 48% of them received antiretroviral therapy (ART), and only 43% of those were viral suppressed; the ART coverage among pregnant women was 94% (see text footnote 1).

Molecular epidemiology and phylogenetic analysis are widely used to characterize HIV-1 epidemics. HIV-1 phylogenetics has been extensively used to characterize the virus transmission networks (Wilkinson et al., 2014; Grabowski et al., 2018) and to reconstruct the viral origin and spread history (Wilkinson et al., 2014). The combination of epidemiological and sociodemographic data with the phylogenetic approaches reveals data on the risk factors associated with the HIV-1 spread (Wolf et al., 2017; Fujimoto et al., 2021), identifies sub-epidemics (Peters et al., 2016), and informs prevention and care interventions (Brenner and Wainberg, 2013).

Molecular–epidemiological and phylogenetic studies of HIV-1 infection in Central Asia, particularly in Kyrgyzstan, are very limited. One of the most recent studies in Kyrgyzstan describes a high rate (range 39–50%) of HIV-1 drug resistance in individuals with treatment failure (Lapovok et al., 2020). The study of HIV-1 sub-subtype A6 [also known as Russian A1 or A-FSU (Foley et al., 2016)] describing the transmission networks in former Soviet Union (FSU) countries concludes that the major driving source of HIV-1 transmission is migrant workers, emphasizing an extensive mobility in the region (80–90% of migrants are from Central Asia) (Aibekova et al., 2018). The study also points at the emerging HIV-1 epidemic among hetero- and homosexual populations, surpassing the parenteral transmission (Aibekova et al., 2018). Drug trafficking from Afghanistan and Tajikistan to southern Kyrgyzstan significantly increased in the mid-1990s, resulting in predominant opiate injections over the homemade drugs (Wolf et al., 2008). Between 1991 and 1999, registered drug use rate in Kyrgyzstan increased sevenfold (Wolf et al., 2008). And in the following 5-year period from 2001 to 2006, Kyrgyzstan had a 15-fold increase of HIV-1 infections, with 76% of cases detected among PWID (Wolf et al., 2008). However, in the last decade, national programs on HIV/AIDS in cooperation with international programs (Thorne et al., 2010; Zhao et al., 2020) have significantly improved HIV-1 care services in the country. While in 2010, only 9% of HIV-infected people received ART, this indicator reached 48% in 2020. In 2010, only 17% HIV-infected people on ART were virally suppressed; this number almost tripled by 2020 (see text footnote 1). By the beginning of 2021, 86.5% of HIV-positive individuals on ART had undetectable viral loads².

In this study, we performed the HIV-1 molecular–epidemiological survey in Kyrgyzstan. To achieve that goal, we analyzed HIV-1 genetic diversity and HIV-1 drug resistance, identified potential transmission clusters, described epidemiological characteristics of studied individuals, and reconstructed the evolutionary history of the virus.

Materials and Methods

Study Population

The blood samples were collected from HIV-positive adults and children who visited local HIV/AIDS clinics of the Ministry of Health of Kyrgyzstan. All the individuals in this study were diagnosed with HIV-1 infection a year prior to the sample collection and were on ART for at least a year. Samples were collected in four provinces (Bishkek, Osh, Jalal-Abad, and Batken) in Kyrgyzstan from 2016 to 2019. Data collected from individuals in Jalal-Abad and Batken were combined (further denoted as JAB) due to the small sample size and geographical closeness of these regions. HIV-1 testing was performed at the study sites according to the national guideline³. Samples were shipped to the Department of Retroviruses, State Research Center of Virology and Biotechnology “Vector” (Koltsovo, Novosibirsk region, Russia) for further testing. Demographic and HIV-related characteristics of individuals were collected at the local healthcare facilities from clinical records. A woman was assigned a “pregnant” status if she was pregnant at the time of her visit for the sample collection. In this study, we refer “children” to individuals 14 years of age and under; “young adults,” between 15 and 24 years old; “adults,” between 25 and 49 years old; and “older adults,” individuals of 50 years of age and above. The study did not recruit individuals from any group; only individuals who visited HIV/AIDS clinics were included. However, samples from children and young adults who had in-hospital acquired HIV-1 infection were particularly collected to investigate those cases.

Amplification of HIV-1 Pol Gene Fragment and Sequencing Analysis

Viral RNA or proviral DNA was extracted using the RealBest DeltaMag kit (Vector-Best, Novosibirsk, Russia) according to the manufacturer’s manual. RNA/DNA was used for the amplification of HIV-1 partial pol gene region coding protease and reverse transcriptase (HXB2 #K03455 reference strain coordinates: 2249–3420). Amplification of the pol gene fragments and sequencing analysis were performed as previously described (Maksimenko et al., 2020).

HIV-1 Subtyping and Drug Resistance Analysis

HIV-1 subtyping was performed using automated HIV-1 subtyping tools REGA v 3.0 (Pineda-Pena et al., 2013), COMET (Struck et al., 2014), and recombinant identification program (RIP) (Siepel et al., 1995). HIV-1 subtypes were also investigated using an approximately maximum-likelihood phylogenetic method using FastTree v2.1.9 (Price et al., 2010) with HIV-1 subtype references from Los Alamos National Laboratory (LANL) HIV Sequence Database, 2020⁴. HIV-1 subtype was assigned if three out of four methods agree. Drug resistance mutations (DRMs) were assessed using Stanford HIV drug resistance database (HIVdb Program) (Shafer, 2006). DRM was considered as major according to Stanford HIV drug resistance database.

Phylogenetic and Cluster Analyses

Phylogenetic analysis was conducted for study and background sequences; background sequences were selected from BLAST as 100 most closely related (with a BLAST score < 1e–50) to each study sequence. Sequences were aligned using MAFFT software (Katoh et al., 2017). Recombination analysis was performed for study and background sequences using RDP4 (Martin et al., 2015). Sites of major DRMs (43 codon positions) were removed. Phylogenetic trees were constructed using IQ-TREE v2 (Minh et al., 2020) under the GTR + G4 + I substitution model. Phylogenetic trees were visualized using interactive Tree of Life (iTOL) (Letunic and Bork, 2021). Monophyletic groups of study sequences with branch support ≥ 80 were considered as a distinct viral lineage (clades). Transmission clusters (≥2 individuals infected by direct/indirect transmission of genetically related HIV-1 variants) were identified by Cluster Picker v1.2.3 (Ragonnet-Cronin et al., 2013) using thresholds of 0.045 of maximum pairwise genetic distance between sequences and a branch support of 90.

Phylodynamic Analysis

Genetic diversity and history reconstruction of the major HIV-1 subtypes were explored using a Bayesian Markov chain Monte Carlo (MCMC) phylogenetic analysis using the BEAST v1.10.4 (Suchard et al., 2018). The temporal structure of the datasets was estimated using TempEst v1.5.3 (Rambaut et al., 2016). Analysis was performed using GTR + G4 + I substitution model, with different combinations of molecular clock models (strict and log-normal uncorrelated relaxed), and coalescent models [constant size, exponential growth, Bayesian Skyline, and Gaussian Markov random field (GMRF) Bayesian Skyride]. The adjustment to the data was estimated using the log marginal likelihood estimation (MLE) using path sampling/stepping-stone sampling (PS/SS). The best-fit model was selected based on the Bayes factor (BF; BF ≥ 3 was considered significant). Two independent MCMC runs were performed for 70 × 10⁶ generations. Convergence of the chains was estimated based on the effective sampling size (ESS; cutoff value over 200 for all the parameters) in Tracer v1.7.1 (Rambaut et al., 2018).

Statistical Analysis

Categorical variables were analyzed using modified Fisher’s test; quantitative variables were analyzed using the Kruskal–Wallis test. The Monte Carlo method (10⁶ simulations) was used for p-value (P) estimation. P correction was performed to control the false discovery rate using the Benjamini–Yekutieli procedure. Statistical analysis was performed in RStudio Team (2020) v1.1.422⁵.

Nucleotide Sequence Accession Numbers

Study sequences were submitted to GenBank under accession numbers MG798935--MG799123, MK228729--MK228833, and MW303524--MW303757 and can be found in the online repository⁶.

Ethics Statement

The study was approved by the Ethical Committee of Research and Manufacturing Association “Prevention Medicine” of Ministry of Health of Kyrgyz Republic (2017). Written consent forms were provided by all the study individuals. For individuals 18 years old and younger, written consent forms were provided by their parents or legal guardians.

Results

Study Population

The study population included 555 individuals who resided in four Kyrgyz provinces (202 individuals from Bishkek, 341 individuals from Osh, and 12 individuals from JAB). Detailed characteristics of the study individuals are presented in Table 1 and Supplementary Table. Men and women were equally represented in the study population (50.8% men vs. 49.2% women), with a median age of 31 years (range 3–72 years old). HIV-1 infection prevailed among individuals in the 25–49 age group. The median time since HIV-1 diagnosis was 8.1 years (range 1–19 years); the median time since HIV-1 diagnosis was longer in Osh than in Bishkek (9 vs. 6.46 years, p < 0.01). Heterosexual transmission was the most prevalent (38.7%), followed by the nosocomial infection (21.4%). In Bishkek, over a half of individuals were infected through heterosexual contacts, and 30% were PWID. In Osh, most of the individuals acquired HIV-1 while getting clinical services. Men were more likely to acquire HIV-1 infection during drug injections (33.5%) and clinical services (25.3%), while women were more likely acquire HIV-1 infection through sexual contacts (72.6%, p < 0.01). When stratified by age group, heterosexual transmission was more prevalent among adults; older adults were more likely to acquire HIV-1 while injecting drugs. All the individuals with nosocomial HIV-1 infection were 24 years old and younger. HIV-1 DRMs to at least two ARV drug classes were more frequent in children and young adults compared with older groups (p < 0.01). Individuals 24 years of age and younger tended to have HIV-1 infection caused by CRF02_AG, while older individuals were more likely to be infected with sub-subtype A6 (p < 0.01). Nine (1.6%) individuals were identified as men who have sex with men (MSM), with a median age of 28 years (range 23–37 years old); seven of those individuals resided in Bishkek province. The median time since HIV-1 diagnosis was 4.6 years (range 4–6 years). Four (44.4%) of nine MSM were infected by minor HIV-1 variants (subtype B, subtype G, URF B/CRF02_AG, and URF B/G); the remaining five MSM had sub-subtype A6 (n = 4) and CRF02_AG (n = 1) HIV-1 infection. None of MSM had HIV-1 DRMs detected.

TABLE 1

Table 1. Demographic and HIV-related characteristics of the 555 study individuals.

HIV-1 Subtyping and Drug Resistance

Genotyping results were successfully obtained for 555 (89.4%) of 621 individuals. The most frequent HIV-1 subtypes were CRF02_AG [n = 332 (59.8%)], CRF63_02A [n = 10 (1.8%)], and sub-subtype A6 [n = 184 (33.2%)]. The remaining 29 (5.2%) sequences represent the minor HIV-1 subtypes and unique recombinant forms (URFs). Detailed HIV-1 subtype distribution is shown in Table 2.

TABLE 2

Table 2. Distribution of HIV-1 subtypes by the study provinces.

Major DRMs were identified in 149 of 555 (26.9%) sequences (Table 2 and Figure 1). Nucleoside reverse-transcriptase inhibitor (NRTI)-resistance mutations were detected in 107 individuals. Non-NRTI (NNRTI)-resistance mutations were detected in 138 individuals. Protease inhibitor (PI)-resistance mutations were detected in four individuals. Resistance mutations to at least two ARV drug classes were identified in 99/149 (66.4%) sequences.

FIGURE 1

Figure 1. Maximum-likelihood phylogenetic trees of HIV-1 pol sequences of CRF02_AG/CRF63_02A (A) and sub-subtype A6 (B). Monophyletic clades (A–D) with the branch support ≥ 80 and containing for over 70% of study sequences are shaded in yellow. JAB, Jalal-Abad and Batken provinces; CRF, circulating recombinant form; MSM, men who have sex with men; PWID, persons who inject drugs; DRM, drug resistance mutation.

Phylogenetic and Cluster Analyses

Phylogenetic trees were constructed for each major HIV-1 subtype separately (Figure 1); sequences of CRF02_AG and CRF63_02A were combined, and a single phylogenetic tree was constructed. Nineteen study sequences (CRF02_AG, n = 11 and A6, n = 8) were excluded from the analysis due to short sequence length (less than 1,000 base pairs). Final datasets included 325 study and 486 background sequences of CRF02_AG/CRF63_02A, and 173 study and 1,056 background sequences of sub-subtype A6. Phylogenetic analyses of CRF02_AG/CRF63_02A and sub-subtype A6 datasets revealed that Kyrgyz sequences intermingle with sequences from the neighboring countries, indicating genetic similarity and the potential common origins with corresponding variants in the Eastern Europe and Central Asia regions. The CRF02_AG/CRF63_02A tree shows the presence of the four well-defined clades composed of the local study sequences; two of those clades (A and B, 97 and 100 branch support values) are predominantly represented by the study sequences from Bishkek (Figure 1A). Two other clades (C and D, branch support values of 84 and 87) primarily include study sequences from Osh. This possibly represents the four distinct locally evolving sub-epidemics of CRF02_AG in Kyrgyzstan. Clades A and B were mainly represented by individuals aged 25–49 years infected by heterosexual and PWID transmission modes. Clades C and D are predominantly presented by children and young adults with the nosocomial HIV-1 infection. No noticeable clades were revealed among sub-subtype A6 study sequences; study samples scattered across the background sequences, indicating multiple independent introduction events of sub-subtype A6 to Kyrgyzstan from the neighboring countries (Figure 1B).

Fifty-seven putative transmission clusters representing 140 individuals (28.1%; 140/498) were detected; 34 clusters were identified among CRF02_AG sequences, 23 clusters among A6 sequences, and one cluster of CRF63_02A sequences (Figure 2). Most of the clusters were pairs (n = 38) and triplets (n = 12). Five clusters of four and one cluster of five sequences were also found. Nine clusters were identified among individuals from different study regions. Major DRMs were detected in 33 (23.6%) clustered individuals; 17 of those individuals had DRMs to at least two ARV dug classes. No statistically significant association between demographic and HIV-related characteristics of the study individuals and phylogenetic clustering was found.

FIGURE 2

Figure 2. Transmission clusters of study HIV-1 pol sequences ofCRF02_AG/CRF63_02A (A) and sub-subtype A6 (B). Clusters were detected in IQ-TREE at maximum 0.045 genetic distances and 90 branch support thresholds. Each figure (circle/square) corresponds to a node on the phylogenetic tree; figure shape, size, and color correspond to an individual’s gender, age, and HIV-1 transmission mode according to the left center footnote. The line drawn between figures indicates that HIV-1 sequences from the respective individuals in phylogenetic tree fell within the clustering thresholds. The black asterisk indicates the HIV-1 drug resistance; the red asterisk indicates the HIV-1 drug resistance to at least two ARV drug classes. Circles with the bold black border correspond to the pregnant women. The solid oval corresponds to a cluster of three individuals infected by CRF63_02A. The dashed lines divide figures corresponding to individuals from the different study regions. JAB, Jalal-Abad and Batken provinces; MSM, men who have sex with men; PWID, persons who inject drugs.

Phylodynamic Analysis

Log-normal uncorrelated relaxed molecular clock model outperforms a strict model based on the BFs in both CRF02_AG/CRF63_02A and sub-subtype A6 datasets. BFs also indicate an advantage of the Bayesian Skyline coalescent model in CRF02_AG/CRF63_02A, and a small advantage of the Gaussian Markov random field (GMRF) Bayesian Skyride coalescent model in the A6 dataset. The estimated median substitution rate for CRF02_AG/CRF63_02A was 2.79 × 10^–3 (95% highest posterior density [HPD]: 1.93 × 10^–3–3.69 × 10^–3), and 2.55 × 10^–3 (95% HPD: 1.69 × 10^–3–3.43 × 10^–3) for sub-subtype A6. The demographic history of HIV-1 CRF02_AG/CRF63_02A and sub-subtype A6 in Kyrgyzstan is presented in Figure 3. CRF02_AG/CRF63_02A Skyline plot reveals a rapid growth phase in ESS between the late 1990s and 2010 followed by a stable phase (Figure 3A). The estimated date of origin of the CRF02_AG epidemic in Kyrgyzstan was January 1997 (95% HPD: February 1986–May 2004). The origin date of CRF63_02A in Kyrgyzstan was estimated as April 2004 (95% HPD: January 1997–April 2009). The origin dates for the four major CRF02_AG clades were estimated as September 2000 (95% HPD: May 1993–August 2005) for clade D, May 2003 (95% HPD: May 1966–March 2008) for clade C, 2005 (95% HPD: March 1999–February 2009) for clade A, and September 2006 (95% HPD: July 2001–May 2010) for clade B. The date of origin for sub-subtype A6 was estimated as June 1995 (95% HPD: August 1985–August 2002). Sub-subtype A6 Skyride plot had a growth phase in ESS in the mid-2010s; the growth rate began to decline somewhat thereafter (Figure 3B).

FIGURE 3

Figure 3. Phylodynamic reconstruction of HIV-1 sCRF02_AG/CRF63_02A (A) and sub-subtype A6 (B) in Kyrgyzstan. (A) The Bayesian Skyline plot was reconstructed for the 325 pol gene sequences of CRF02_AG and 10 sequences of CRF63_02A. The horizontal bold line indicates effective population size through time; blue-shaded area represents the 95% highest posterior density (HPD). The vertical bold black line indicates the estimated origin date for CRF02_AG [January 1997 (95% HPD: February 1986–May 2004)]; the vertical dashed line indicates the estimated origin date for CRF63_02A [April 2004 (95% HPD: January 1997–April 2009)]. (B) The Gaussian Markov random field (GMRF) Bayesian Skyride plot was reconstructed for the 173 pol gene sequences of sub-subtype A6. The horizontal bold line indicates effective population size through time; pink-shaded area represents the 95% HPD. The vertical bold black line indicates the estimated origin date for sub-subtype A6 [June 1995 (95% HPD: August 1985–August 2002)].

Discussion

In this study, we performed molecular–epidemiological analysis of HIV-1 in Kyrgyzstan in four provinces, two of which (Bishkek and Osh) are most severely affected by HIV-1. Our results showed the predominant circulation of two major HIV-1 variants—CRF02_AG and sub-subtype A6. A third of HIV-infected individuals in the study had major HIV-1 DRMs: two-third of those had DRMs to at least two ARV drug classes. Phylogenetic analysis revealed a well-defined structure of CRF02_AG indicating the locally evolving sub-epidemics. The lack of a well-defined phylogenetic structure was observed for the A6 sub-subtype. Small distinct HIV-1 transmission clusters were detected among study sequences. Clustering was not associated with any individual characteristics. DRMs were detected in 26.9% of clustered individuals.

HIV-1 epidemic in Kyrgyzstan is characterized by an overall low HIV-1 infection prevalence rate, a high HIV-1 prevalence rate among key populations (Deryabina et al., 2019), and a high proportion of HIV-positive children (Mansfeld et al., 2015). Our study found that most of the HIV-positive older adults were PWID; heterosexual transmission was more common in the 25–49 age group. This may indicate the shift of the HIV-1 epidemic to the general population from key groups such as PWID or change of substance consumption to non-injection drugs. Our study also found that over a half of children and young adults in Osh acquired healthcare-associated HIV-1 infection. Earlier studies described several outbreaks of hospital-acquired HIV-1 among children in Osh province in 2007 (Thorne et al., 2010) and 2011–2012 (Mansfeld et al., 2015), with the total number of over 300 reported cases. Although healthcare-associated HIV-1 infections in high-income countries are extremely low, risk of nosocomial HIV-1 infections in the resource-limited countries remains high (Ganczak and Barss, 2008; Myburgh et al., 2020). We also identified significant differences in the transmission modes between genders. Women were more likely to be infected thought the sexual contacts; PWID was a dominant transmission mode among men, followed by nosocomial infection. Similar epidemiological characteristics of HIV-positive individuals were described in the previous report (Zhao et al., 2020).

The first HIV-1 case in Kyrgyzstan was identified in 1987 in a foreigner [34]. The first local HIV-1 case was registered in 1996 (Roth et al., 2012), and HIV-1 cases in Kyrgyzstan are steadily increasing since then, like in other countries of Central Asia and Eastern Europe. CRF02_AG predominantly circulates in West Africa (Bbosa et al., 2019), but since 1999, this recombinant form is constantly detected in Central Asian countries (Laga et al., 2015; Aibekova et al., 2018). CRF02_AG circulating in Central Asia is genetically distant from that from African countries and phylogenetically represents well-supported distant clade (Mir et al., 2016). Our study showed the presence of the four main CRF02_AG region-specific lineages circulating in Kyrgyzstan. Phylodynamic reconstruction of HIV-1 epidemics in Kyrgyzstan revealed the date of origin of CRF02_AG as January 1997. Osh-specific clades increased earlier than did the Bishkek-specific clades. Analysis of the characteristics of the study individuals also showed that people from Osh had HIV-1 infection for a longer period than those from Bishkek. The CRF02_AG rise in Osh coincided with the active drug trafficking from Afghanistan through the so-called ‘‘North route.’’ The main traffic routes go to the South Kyrgyzstan (Osh province) through the Tajikistan and Uzbekistan and then spread across Kyrgyzstan and further to Kazakhstan, Russian, and other European countries⁷. Most likely, this viral variant was introduced to Kyrgyzstan by drug trafficking [Mir et al., 2016]. Bishkek-specific clades were predominantly represented by adults who acquired HIV-1 infection through heterosexual contacts and while injection of drugs, indicating an extensive mixing between high-risk and heterosexual populations. Osh-specific clades are represented by children and young adults who acquired healthcare-associated HIV-1 infection and represented by in-hospital HIV-1 outbreak sequences. These age groups play an important role in the local HIV-1 sub-epidemic, and targeted HIV-1 care interventions should be implemented to limit HIV-1 spread and improve healthcare services. HIV-1 sub-subtype A6 epidemic in FSU countries is characterized by monophyletic phylogenetic tree structure, suggesting a common ancestor for all the viruses (Díez-Fuertes et al., 2015). Our results revealed a significant geographic dissemination of the sub-subtype A6 across the FSU countries with the lack of country-specific clades. HIV-1 sub-subtype A6 phylogenetics discovered multiple introductions of the virus to Kyrgyzstan, which are most likely associated with the intense migration between Kyrgyzstan and the neighboring countries. The origin date for sub-subtype A6 in Kyrgyzstan was estimated as June 1995. The HIV-1 effective population size of these subtypes had an initial growing phase with the further stable phase for CRF02_AG/CRF63_02A and decline phase for sub-subtype A6. CRF63_02A (formally known as CRF02_AG recombinant) was initially described in 2012 as a viral variant that caused an outbreak in the Asian part of Russia in 2006 (Baryshev et al., 2012). In our study, a very limited number of CRF63_02A sequences were detected. A small number (1.6%) of HIV-positive MSM was found in the study. Almost a half of those individuals were infected by minor HIV-1 variants (subtypes B and G, and URFs B/CRF02_AG and B/G); the remaining five MSM had sub-subtype A6 or CRF02_AG HIV-1 infection. Despite the very limited number of MSM in our study, these data indicate the mixing of HIV-1 transmissions between MSM and general population. Also, the high HIV-1 genetic diversity among MSM could potentially be one of the driving forces for further increase of HIV-1 variety in Kyrgyzstan. Data of HIV-1 prevalence among MSM in Central Asia, and particularly in Kyrgyzstan, are very limited due to criminalization of same-sex relationships (Thorne et al., 2010; Wirtz et al., 2013). However, previous studies reported HIV-1 prevalence rate among MSM at around 1% (Thorne et al., 2010; El-Bassel et al., 2013). Further studies are needed to better understand HIV-1 epidemiology and HIV-related risk behavior among this high-risk group in Kyrgyzstan. This study showed that groups of at-risk people could contribute to the HIV-1 epidemic progression by maintaining new infections within the risk group or by linking infections between different populations.

Our study has several limitations. Children and young adults were sampled more intensely than older individuals, which could have biased some of the observed numbers. In this study, samples from children and young adults were deliberately collected from those who potentially acquired HIV-1 infection during in-hospital outbreaks. Also, this study provides information on HIV-positive pregnant women in Bishkek province only. Data on pregnancy status from other provinces were not available. In this study, we detected a small number of transmission clusters sized two to four individuals due to a small sampling fraction. The study dataset represents only around 6% (555 HIV-positive study individuals of an estimated 9,200 HIV-infections in 2020 (see text footnote 1) of the total HIV-1 detected cases in Kyrgyzstan. Major DRMs were detected in 23.6% of individuals in the HIV-1 transmission clusters. The study as performed does not power to reveal whether DRMs occurred due to the therapy failure or whether originating from the transmission of the drug-resistant virus. HIV-positive individuals not linked to care were not included in the study, since only individuals who visited healthcare facilities were included.

To our knowledge, this is the first study that provides comprehensive data on the HIV-1 epidemic in Kyrgyzstan. We identified complex patterns of HIV-1 phylogenetics among the key and general populations and observed high levels of HIV-1 genetic diversity and drug resistance. Further studies are needed, accompanied by extensive public health interventions to limit HIV-1 spread and improve future HIV-1 care services.

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.

Ethics Statement

The study was approved by Ethical Committee of Research and Manufacturing Association “Prevention Medicine” of Ministry of Health of Kyrgyz Republic. Written consent forms were provided by all the study individuals. For individuals 18 years old and younger, written consent forms were provided by their parent or legal guardian.

Author Contributions

MS, AT, and NG planned and designed the study. DZ, IO, TN, MG, and VI performed HIV genotyping and collected the sequences. AT submitted the sequences. MS performed the analysis of the epidemiological data, performed the phylogenetic and phylodynamic analyses, produced the illustrations, and wrote the manuscript. NG supervised the project and edited the manuscript. IM performed the statistical analysis. UC, EN, UM, ZA, NA, AB, and UK collected the epidemiological data and assisted with sample collection at the study sites. All authors participated in the final review of the manuscript.

Funding

This study was supported by the Directives of the Russian Government Nos. 2314-p dated November 14, 2015 and No. 2656-p dated December 1, 2018.

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.

Acknowledgments

We thank the study individuals for their contributions to the study. We also thank the staff at the study HIV/AIDS clinics for sample management and laboratory testing.

Supplementary Material

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

Footnotes

1. ^ https://www.unaids.org/en/regionscountries/countries/kyrgyzstan
2. ^ http://aidscenter.kg/statistika/
3. ^ http://med.kg/images/MyFiles/KP/2018/sbornik_KP_VISH_prikaz_903_10102017.pdf
4. ^ http://www.hiv.lanl.gov/
5. ^ http://www.rstudio.com/
6. ^ https://www.ncbi.nlm.nih.gov/genbank/
7. ^ https://www.unodc.org/documents/publications/NR_Report_21.06.18_low.pdf

References

Aibekova, L., Foley, B., Hortelano, G., Raees, M., Abdraimov, S., Toichuev, R., et al. (2018). Molecular epidemiology of HIV-1 subtype A in former Soviet Union countries. PLoS One 13:e0191891. doi: 10.1371/journal.pone.0191891

PubMed Abstract | CrossRef Full Text | Google Scholar

Baryshev, P. B., Bogachev, V. V., and Gashnikova, N. M. (2012). Genetic characterization of an isolate of HIV type 1 AG recombinant form circulating in Siberia, Russia. Arch. Virol. 15, 2335–2341. doi: 10.1007/s00705-012-1442-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bbosa, N., Kaleebu, P., and Ssemwanga, D. (2019). HIV subtype diversity worldwide. Curr. Opin. HIV AIDS 14, 153–160. doi: 10.1097/COH.0000000000000534

PubMed Abstract | CrossRef Full Text | Google Scholar

Brenner, B. G., and Wainberg, M. A. (2013). Future of phylogeny in HIV prevention. J. Acquir. Immune Defic. Syndr. 63 Suppl 2, 248–254. doi: 10.1097/QAI.0b013e3182986f96

PubMed Abstract | CrossRef Full Text | Google Scholar

Deryabina, A. P., Patnaik, P., and El-Sadr, W. M. (2019). Underreported injection drug use and its potential contribution to reported increase in sexual transmission of HIV in Kazakhstan and Kyrgyzstan. Harm Reduct. J. 16:1. doi: 10.1186/s12954-018-0274-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Díez-Fuertes, F., Cabello, M., and Thomson, M. M. (2015). Bayesian phylogeographic analyses clarify the origin of the HIV-1 subtype A variant circulating in former Soviet Union’s countries. Infect. Genet. Evol. 33, 197–205. doi: 10.1016/j.meegid.2015.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Bassel, N., Strathdee, S. A., and El Sadr, W. M. (2013). HIV and people who use drugs in central Asia: confronting the perfect storm. Drug Alcohol Depend 132 Suppl 1, S2–S6. doi: 10.1016/j.drugalcdep.2013.07.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Foley, B. T., Leitner, T., Paraskevis, D., and Peeters, M. (2016). Primate immunodeficiency virus classification and nomenclature: review. Infect. Genet. Evol. 46, 150–158. doi: 10.1016/j.meegid.2016.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimoto, K., Bahl, J., Wertheim, J. O., Del Vecchio, N., Hicks, J. T., Damodaran, L., et al. (2021). Methodological synthesis of Bayesian phylodynamics, HIV-TRACE, and GEE: HIV-1 transmission epidemiology in a racially/ethnically diverse Southern U.S. context. Sci. Rep. 11:3325. doi: 10.1038/s41598-021-82673-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganczak, M., and Barss, P. (2008). Nosocomial HIV infection: epidemiology and prevention – a global perspective. AIDS Rev. 10, 47–61.

Google Scholar

Grabowski, M. K., Herbeck, J. T., and Poon, A. F. Y. (2018). Genetic cluster analysis for HIV prevention. Curr. HIV/AIDS Rep. 15, 182–189. doi: 10.1007/s11904-018-0384-1

PubMed Abstract | CrossRef Full Text | Google Scholar

HIV Sequence Database (2020). HIV Sequence Database. Available online at: http://www.hiv.lanl.gov/ (accessed November 15, 2020).

Google Scholar

Katoh, K., Rozewicki, J., and Yamada, K. D. (2017). MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 20, 1160–1166. doi: 10.1093/bib/bbx108

PubMed Abstract | CrossRef Full Text | Google Scholar

Laga, V., Lapovok, I., Kazennova, E., Ismailova, A., Beisheeva, N., Asybalieva, N., et al. (2015). The genetic variability of HIV-1 in Kyrgyzstan: the spread of CRF02_AG and subtype A1 recombinants. J. HIV AIDS 1. doi: 10.16966/2380-5536.106

CrossRef Full Text | Google Scholar

Lapovok, I. A., Saleeva, D. V., Lopatukhin, A. E., Kirichenko, A. A., Murzakova, A. V., Bekbolotov, A. A., et al. (2020). HIV-1 genetic characteristics in patients having experienced virological failure of therapy in the Kyrgyz Republic in 2017–2018. Epidemiol. Infect. Dis. 3, 105–111. doi: 10.18565/epidem.2020.10.3.105-11

CrossRef Full Text | Google Scholar

Letunic, I., and Bork, P. (2021). Interactive tree of life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 49, W293–W296. doi: 10.1093/nar/gkab301

PubMed Abstract | CrossRef Full Text | Google Scholar

Maksimenko, L. V., Totmenin, A. V., Gashnikova, M. P., Astakhova, E. M., Shudarnov, S. E., Ostapova, T. S., et al. (2020). Genetic diversity of HIV-1 in Krasnoyarsk Krai: area with high levels of HIV-1 recombination in Russia. Biomed. Res. Int. 10:9057541. doi: 10.1155/2020/9057541

PubMed Abstract | CrossRef Full Text | Google Scholar

Mansfeld, M., Ristola, M., and Likatavicius, G?. (2015). HIV Programme Review in Kyrgyzstan. Copenhagen: WHO Regional Office for Europe.

Google Scholar

Martin, D. P., Murrell, B., Golden, M., Khoosal, A., and Muhire, B. (2015). RDP4: detection and analysis of recombination patterns in virus genomes. Virus Evol. 1:vev003. doi: 10.1093/ve/vev003

PubMed Abstract | CrossRef Full Text | Google Scholar

Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., Woodhams, M. D., von Haeseler, A., et al. (2020). IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534. doi: 10.1093/molbev/msaa015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ministry of Health of Kyrgyz Republic (2017). Clinical Protocols on HIV Infection for Outpatient and Inpatient Services. Bishkek. Available online at: http://med.kg/images/MyFiles/KP/2018/sbornik_KP_VISH_prikaz_903_1010 2017.pdf (accessed May 20, 2021).

Google Scholar

Mir, D., Jung, M., Delatorre, E., Vidal, N., Peeters, M., and Bello, G. (2016). Phylodynamics of the major HIV-1 CRF02_AG African lineages and its global dissemination. Infect. Genet. Evol. 46, 190–199. doi: 10.1016/j.meegid.2016.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Myburgh, D., Rabie, H., Slogrove, A. L., Edson, C., Cotton, M. F., and Dramowski, A. (2020). Horizontal HIV transmission to children of HIV-uninfected mothers: a case series and review of the global literature. Int. J. Infect. Dis. 98, 315–320. doi: 10.1016/j.ijid.2020.06.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, P. J., Pontones, P., Hoover, K. W., Patel, M. R., Galang, R. R., Shields, J., et al. (2016). HIV infection linked to injection use of oxymorphone in Indiana, 2014–2015. N. Engl. J. Med. 375, 229–239. doi: 10.1056/NEJMoa1515195

PubMed Abstract | CrossRef Full Text | Google Scholar

Pineda-Pena, A. C., Faria, N. R., Imbrechts, S., Libin, P., Abecasis, A. B., Deforche, K., et al. (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

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, M. N., Dehal, P. S., and Arkin, A. P. (2010). FastTree 2 – approximately maximum-likelihood trees for large alignments. PLoS One 5:e9490. doi: 10.1371/journal.pone.0009490

PubMed Abstract | CrossRef Full Text | Google Scholar

Ragonnet-Cronin, M., Hodcroft, E., Hué, S., Fearnhill, E., Delpech, V., Leigh Brown, A. J., et al. (2013). Automated analysis of phylogenetic clusters. BMC Bioinformatics 6:317. doi: 10.1186/1471-2105-14-317

PubMed Abstract | CrossRef Full Text | Google Scholar

Rambaut, A., Drummond, A. J., Xie, D., Baele, G., and Suchard, M. A. (2018). Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904. doi: 10.1093/sysbio/syy032

PubMed Abstract | CrossRef Full Text | Google Scholar

Rambaut, A., Lam, T. T., Max Carvalho, L., and 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, S., Elfving, R., and Khan, A. R. (2012). HIV/AIDS Vulnerabilities in Regional Transport Corridors in the Kyrgyz Republic and Tajikistan. Manila: Asian Development Bank.

Google Scholar

RStudio Team (2020). RStudio: Integrated Development for R. RStudio. Boston, MA: RStudio Team.

Google Scholar

Shafer, R. W. (2006). Rationale and uses of a public HIV drug-resistance database. J. Infect. Dis. 194(Suppl. 1) 51–58. doi: 10.1086/505356

PubMed Abstract | CrossRef Full Text | Google Scholar

Siepel, A. C., Halpern, A. L., Macken, C., and Korber, B. T. (1995). A computer program designed to screen rapidly for HIV type 1 intersubtype recombinant sequences. AIDS Res. Hum. Retroviruses 11, 1413–1416. doi: 10.1089/aid.1995.11.1413

PubMed Abstract | CrossRef Full Text | Google Scholar

Struck, D., Lawyer, G., Ternes, A. M., Schmit, J. C., and Bercoff, D. P. (2014). COMET: adaptive context-390based modeling for ultrafast HIV-1 subtype identification. Nucleic Acids Res. 39142:e144. doi: 10.1093/nar/gku739

PubMed Abstract | CrossRef Full Text | Google Scholar

Suchard, M. A., Lemey, P., Baele, G., Ayres, D. L., Drummond, A. J., and Rambaut, A. (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4:vey016. doi: 10.1093/ve/vey016

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorne, C., Ferencic, N., Malyuta, R., Mimica, J., and Nieniec, T. (2010). Central Asia: hotspot in the worldwide HIV epidemic. Lancet Infect. Dis. 10, 479–488. doi: 10.1016/S1473-3099(10)70118-3

CrossRef Full Text | Google Scholar

Wilkinson, E., Engelbrecht, S., and de Oliveira, T. (2014). Detection of transmission clusters of HIV-1 subtype C over a 21-year period in Cape Town, South Africa. PLoS One 9:e109296. doi: 10.1371/journal.pone.0109296

PubMed Abstract | CrossRef Full Text | Google Scholar

Wirtz, A. L., Kirey, A., Peryskina, A., Houdart, F., and Beyrer, C. (2013). Uncovering the epidemic of HIV among men who have sex with men in Central Asia. Drug Alcohol Depend 132(Suppl. 1) S17–S24. doi: 10.1016/j.drugalcdep.2013.06.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, E., Herbeck, J. T., Van Rompaey, S., Kitahata, M., Thomas, K., Pepper, G., et al. (2017). Phylogenetic evidence of HIV-1 transmission between adult and adolescent men who have sex with men. AIDS Res. Hum. Retroviruses 33, 318–322. doi: 10.1089/AID.2016.0061

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, D., Elovich, R., Boltaev, A., and Pulatov, D. (2008). HIV in Central Asia: Tajikistan, Uzbekistan and Kyrgyzstan. Public Health Aspects of HIV/AIDS in Low- and Middle-Income Countries. New York, NY: Springer.

Google Scholar

Zhao, F., Benedikt, C., and Wilson, D. (2020). Tackling the World’s Fastest-Growing HIV Epidemic: More Efficient HIV Responses in Eastern Europe and Central Asia. Human Development Perspectives. Washington, DC: World Bank.

Google Scholar