Following China’s open-door policy, the use of illicit substances emerged, which was associated with global drug-trafficking activities in the 1980s^[1]. The proximity of the “Golden Triangle” (the border areas of Thailand, Myanmar, and Laos) facilitates drug trafficking from the origin of the drug to Yunnan and Guangxi, two provinces in southwest China, or from Vietnam to Guangxi, and then to other neighboring provinces^[2,3]. In 1996, the first HIV case in Guangxi was identified among people who inject drugs (PWID)^[4]. Since then, Guangxi has become an HIV-1 hotspot, with many HIV-infected PWID^[5]. Until 2006, IDUs were a major contributor to the spread of HIV in China, although sexual transmission was subsequently identified as the main route of HIV transmission^[6]. From 2012 to 2019, 28.0% of newly diagnosed cases of HIV infected via IDU had baseline CD4 cell counts greater than or equal to 500 cells/μL, indicating that these individuals may have been infected within the past 2 years^[7]. Since 2010, the HIV prevalence among new entrants tested for HIV at detoxification centers has been increasing^[8], which means that HIV transmission via IDU still exists^[7].

Previous evidence has said that HIV is characterized by mutation and recombination, which leads to high genetic diversity and complex heterogeneity^[9]. HIV-1 subtype B’ (Tai B) strains and C were co-circulated among IDU in Yunnan Province in the early 1990s^[10-12]. Subsequently, CRF07_BC, CRF08_BC, and a large number of unique recombinant forms began to circulate widely among IDUs in China^[13,14]. CRF01_AE was found in populations at risk of sexual transmission and IDUs in Yunnan and Guangxi^[15-17], bordering Thailand, Myanmar, and Vietnam, and quickly became the most widespread HIV-1 strain in the national molecular epidemiological survey^[6]. Guangxi’s HIV epidemic had become one of the worst in China as the initial epidemic was driven by IDUs in the 1990s. In Guangxi, HIV subtypes B′ (Tai B), C, D, and CRF01_AE were identified in 1996^[17]. HIV subtype B/C recombinants were found among IDUs in Yunnan Province and Baise City, Guangxi^[5]. During this period, the CRF01_AE strain was found among IDUs in Pingxiang and Nanning cities in Guangxi, which are located along the main drug-trafficking route from Vietnam to China^[18]. CRF07_BC and CRF08_BC, originally formed in Yunnan Province, were first detected among IDUs in Guangxi in 1997^[19]. In contrast, the first CRF01_AE strain identified among IDUs in Guangxi had a significant clustering with strains found in Northern Vietnam^[20]. Recent studies have found that HIV dominant subtypes among IDUs in China are CRF07-BC, followed by CRF01-AE and CRF08-BC^[21].

In Guangxi, HIV-1 transmission via IDU has not yet been eliminated, and the epidemic is characterized by substantial viral diversity and distinct molecular properties. However, there have been few studies on the regional HIV epidemic among PWIDs in recent years. A detailed and comprehensive analysis of the geographic distribution of viral genotypes is needed. In this study, to determine the source of infections, we provide evidence of the genotypes and the temporal-geographical origins of HIV transmission among IDUs in Guangxi. Furthermore, we spatialized the HIV-1 genetic transmission network, with an explanation of spatial distribution and spatial clustering patterns of phylogenetic clustering rates and spatial links of HIV-1 among cities, which helps to better understand the dynamic characteristics toward the IDU-related transmission of HIV-1 subtypes at the regional level and may contribute to developing strategies to reduce the risk of HIV transmission via IDU.

The blood samples from 535 confirmed HIV-1-positive IDUs were continuously collected between January 2009 and December 2019 through voluntary counseling and testing sites in Guangxi, China. We performed a BLAST search using the sequences and selected the top five sequences with the highest homology from publicly available HIV databases at Los Alamos National Laboratory^[22]. Seven hundred fifty reference sequences from the HIV database were selected for analysis after excluding repeated sequences (accessed May 2021). Overall, 1,285 sequences were included for subsequent analyses.

Sequencher 5.1 software was used to edit all sequences, and BioEdit 7.1 software was used to align them. The reference sequences were obtained from the Los Alamos National Laboratory database and covered the major HIV-1 subtypes and circulating recombinant forms (CRFs). Phylogenetic tree analyses were performed using the neighbor-joining method based on the Kimura2-parameter model with 1,000 bootstrap replicates by using MEGA7.0 software^[23].

One thousand two hundred eighty-five aligned sequences, including reference and subjects’ sequences, were entered into the HyPhy software to calculate genetic distances, and Tamura-Nei93 pairwise genetic distances were calculated for all pairs of sequences. Two sequences with a genetic distance of 0.6% were identified as potential transmission partners. Cytoscape V3.5.1 was used to process the network data and construct the transmission network. As previously mentioned, the characteristics of the network, including nodes (individuals in the network), edges (the link between two nodes, representing the potential transmission relationship between the two individuals), degrees (the number of edges linking one node to other nodes), network sizes (the number of individuals in a cluster), and clusters (groups of linked sequences), were determined^[24]. The Cytoscape V3.5.1 software was developed by Institute for Systems Biology. The network data were visualized using ArcGIS Desktop 10.7, Version: 10.7.0.10450 (Environmental Systems Research Institute)^[25].

For accurate analysis and better display of phylogenetic trees, 516 pol sequences were divided into three data sets according to genotypes CRF01_AE, CRF07_BC, and CRF08_BC, with 152, 57, and 307 sequences, respectively. Since only one sequence belonged to subtype C and CRF57_BC, no phylogenetic tree was required to determine the transmission source. Furthermore, additional 666 unique pol sequences identified by BLAST as CRF01_AE, CRF07_BC, and CRF08_BC, were divided into data sets with 368, 134, and 164 sequences, respectively. Thus, 1,182 (516 + 666) pol sequences were used to determine the sources of HIV-1 transmission strains.

To reconstruct the spatial dynamics and estimate the sources of the strains, a Bayesian discrete phylogeographic approach was performed using Markov chain Monte Carlo (MCMC) runs of 300 million generations with BEAST v.1.8.4 under a Bayesian skygrid demographic model. The first 10%–30% of the states from each run were discarded as the burn-in^[26,27]. All three data sets were analyzed using a general time-reversible model specifying a gamma distribution based on a neighbor-joining tree as a prior on each relative substitution rate and a relaxed uncorrelated lognormal molecular clock model to infer the timescale of HIV evolution with a gamma distribution prior on the mean clock rate (shape = 0.001, scale = 1,000)^[28]. The Bayesian MCMC output was analyzed using Tracer v1.6 (http://beast.bio.ed.ac.uk/Tracer). Maximum clade credibility (MCC) trees were generated in TreeAnnotator v1.8.4. Finally, the most probable temporal, geographical origin, and risk category of HIV transmission strains for the study participants were inferred according to the output of the posterior of Bayesian estimation and visualized on MCC trees using the FigTree software v1.4.3 (http://beast.bio.ed.ac.uk). The tree figure editing was conducted using Adobe Illustrator 2020 (V 24.0.1.341).

Sequences with a posterior probability over 0.7 were defined as local circulating strains and risk categories in Guangxi, and the posterior probability refers to the support of locations at the nodes. The provinces of origin of HIV transmission strains in our study samples were then tabulated according to different HIV subtypes.

The edges between the two nodes were calculated to explore factors associated with potential transmission. Comparisons between individuals were based on the following factors: year of diagnosis, residence region, age, marital status, educational level, ethnicity, HIV subtype, and time of infection. Significant differences in categorical variables were analyzed using chi-square or Fisher’s exact tests. P < 0.05 was considered statistically significant.

Table 1 provides a summary of the 535 HIV-positive IDUs samples. The average age of these individuals was 42.53 ± 8.13 years, of which 66.7% (357/535) were over 35 years old. There were 489 (91.4%) women and 46 (8.6%) men. Of these IDUs, 38.9% (208/535) were permanent residents in Nanning, 33.6% (327/535) were from Qinzhou, and 27.5% were from other cities in Guangxi. The ethnicities of the participants were mainly Han (411, 76.8%) and Zhuang (105, 19.6%), and the education level was mainly junior school or below (484, 90.5%). In terms of marital status, 61.3% (328/535) were single, and 35.7% (191/535) were married. For CD4⁺ T-cell count, 110 (20.6%) individuals had a CD4⁺ T-cell count ≥ 500 cells/μL, 163 (30.5%) had a CD4⁺ T-cell count < 250 cells/μL, and 194 (36.3%) between 250 and 500 cells/μL. Of all participants, 320 subjects (59.8%) were ART naïve (Table 1).

As shown in Figure 1, 535 samples were sequenced and genotyped. HIV-1 subtypes were mainly CRF08_BC, CRF01_AE, and CRF07_BC. A few other subtypes/CRFs were C and CRF57_BC. The 152 CRF01_AE sequences were divided into four clusters: CRF01_AE Cluster 1 (82/152, 53.9%), CRF01_AE Cluster 2 (33/152, 21.7%), CRF01_AE Cluster 3 (3/152, 2.0%), and CRF01_AE Cluster 4 (2/152, 1.3%) (Figure 1, Table 2). Seventeen (3.2%) URFs were discovered, mainly between CRF_01AE and CRF_07/08_BC (Supplementary Figure S1, Supplementary Table S1, available in www.besjournal.com). There were significant differences in the distribution of genotypes by the time of HIV diagnosis, participants’ region, age, education, and treatment (Table 1).

Figure 1.  Maximum-likelihood trees of 535 HIV-1 pol sequences sampled from IDUs in Guangxi, China. The stability of each node was assessed by bootstrap analyses with 1,000 replicates. The red spots and green spots indicate the sequences of this study and references from the HIV database, respectively.

At a genetic distance threshold of 0.6%, transmission network analysis found 34.3% (441/1,285) genetically linked individuals, forming 128 clusters and ranging in size from 2 to 80. The largest cluster consisting of 595 links was comprised predominantly of subtype CRF08_BC, which clustered more easily than other subtypes, indicating that CRF08_BC tended to group into a large cluster. In the network, 121 subjects (27.4%) were from Guangxi, and 320 (72.6%) were from other provinces (Figure 2). To explore factors associated with the potential transmission, 535 Guangxi HIV-positive IDUs samples were included in the analysis. The chi-square test revealed significant differences in the transmission network between subtypes (P = 0.005), age (P < 0.001), region (P < 0.001), year of diagnosis (P < 0.001), and treatment status (P < 0.001) (Table 3).

Figure 2.  Region and transmission route associated with HIV-1 transmission clusters. The network was constructed using R script with the network package in R version 4.0.2 software. Shape codes the reported transmission route, and colors denote regions for each individual. HET, heterosexual; IDU, intravenous drug user; MSM, men who have sex with men.

To explore the geographic dimension of HIV-1 transmission among IDUs, we constructed transmission routes using the R language. Guangxi was linked with Vietnam (2.7%) and 11 provinces in China, including Guangdong (43.5%), Yunnan (16.7%), and others (37.1%). There were 99 (22.3%) intercity links from Guangxi to other provinces, with Guangdong (44, 9.9%) having the strongest link. In addition, 248 (55.8%) intercity links were from other provinces to Guangxi, of which Guangdong (128, 28.8%) had the strongest link, followed by Yunnan (69, 15.5%) (Figure 3).

Figure 3.  The geographic dimension of HIV-1 transmission. HIV-1 transmission between Guangxi and other provinces. The line width indicates the number of individuals.

Considering that CRF08_BC, CRF01_AE, and CRF07_BC were the dominant strains, other subtypes were excluded from the analysis due to small sample sizes (n < 5).

The sources of HIV-1 transmission strains in recently infected individuals in Guangxi were estimated by the reconstruction of Bayesian discrete phylogeographic approaches under a Bayesian skygrid demographic model. Based on the criteria of posterior probability over 0.7, the transmission source of 516 samples could be determined. As shown in Table 4, 483 (93.6%) samples were derived from strains already circulating in Guangxi, with the remaining 33 (6.4%) derived from strains circulating outside Guangxi. HIV transmission sources varied by subtypes and education levels. Subtype CRF01_AE had the most diversified source (Table 5).

The introduction time for CRF01_AE strains in Guangxi IDUs with effective population size circulating may be 1994.25 (1992.97–1995.58), from Vietnam but originated in Thailand. Two distinct phylogenetic clusters of CRF01_AE were identified. Cluster strains were placed within a large phylogenetic group consisting of CRF01_AE strains that originated in Thailand and then spread to Guangxi through Vietnam. This cluster contained CRF01_AE strains circulating among IDUs in Guangxi, Guangdong, and Yunnan. Another cluster originated in Thailand and then spread to Guangxi through Guangdong and formed three clusters of Guangxi local circulating strains. CRF07_BC strains probably came from Guangdong in 1995.23 (1993.64–1996.64). CRF07_BC strains in Guangxi may originate in Yunnan and then spread to Guangxi through Guangdong. It is worth noting that some individuals may be associated with the epidemic among MSMs in Guangdong. CRF08_BC strains, mainly directly originating from Yunnan, probably spread to Guangxi in 1992.29 (1989.05–1995.36), forming large local clusters among IDUs (Supplementary Figure S2, available in www.besjournal.com).

Figure S2.  This figure illustrated time-scaled phylogeographic history of three HIV-1 subtypes (CRF01_AE, CRF08_BC, CRF07_BC) for IDU in Guangxi. Background colours around branches represented the most probable origin province of the parental node of each branch. The names of the origin provinces were indicated by the colour of the legend in the upper left corner. The values of probability inferred by Bayesian skygrid demographic model estimation for the most probable ancestral location were specified at the corresponding nodes. The probable time of introduction for these three key subtypes with an effective population size circulating in Guangxi IDU was indicated at key nodes. Red and blue dots at the ends of branches denoted the sequences of study participants and references, respectively. The Bayesian MCC trees were visualised through the FigTree software v1.4.3 (http://beast.bio.ed.ac.uk). The figure editing was using Adobe Illustrator 2020 (V 24.0.1.341).

HIV transmission risk categories for 516 sequences with CRF08_BC, CRF01_AE, and CRF07_BC were visualized on the MCC tree based on a posterior probability criterion exceeding 0.7 (Supplementary Figure S2). Of these, 92.1% of sequences were more frequently linked to others reporting IDU risk. The percentage of CRF08_BC sequences linked to IDU risk was as high as 97.7%, followed by CRF07_BC (94.7%) and CRF01_AE (80.1%). In addition, some were also linked to reporting heterosexual (7.2%) and MSM (0.6%) risks. Compared with the CRF08_BC and CRF07_BC groups, the CRF01_AE group had a higher proportion of IDUs (18.5%) linked to those reporting heterosexual risk (Supplementary Figure S3, available in www.besjournal.com).

Figure S3.  This figure illustrated HIV transmission risk category link of three HIV-1 subtypes (CRF01_AE, CRF07_BC, CRF08_BC) for IDU in Guangxi. The transmission risk category of the strains is color-coded as shown in the inset. Red and blue dots at the ends of branches denoted sequences of study participants and references, respectively. Background colors around branches represent the most probable origin of the parental node of each branch. The Bayesian MCC trees were visualized through the FigTree software v1.4.3 (http://beast.bio.ed.ac.uk). The figure editing was using Adobe Illustrator 2020 (V 24.0.1.341).

This study investigated the source and genetic transmission patterns of HIV-1 strains circulating among IDUs in Guangxi using molecular epidemiological approaches. Although a previous study discovered that HIV-1 transmission in Guangxi predominantly occurs from strains already circulating locally, the small sample size in a previous study failed to indicate the current HIV/AIDS status of IDUs in Guangxi^[23]. In recent years, HIV infection due to drug use still disproportionately affects the AIDS epidemic, and the transmission patterns of this high-risk group still need to be researched. To the best of our knowledge, this is the first large-scale study of HIV genetic transmission among drug users in Guangxi using molecular epidemiological methods. We found that HIV-1 transmission among recently diagnosed IDUs was predominantly derived from locally circulating strains and formed epidemic clusters in Guangxi. The major HIV-1 subtype strains among IDUs in Guangxi were CRF08_BC, CRF01_AE, and CRF07_BC, which came from different regions, respectively, rather than a single source. The most probable source of CRF01_AE strains was either directly from Vietnam to Guangxi in 1994.25 (1992.97–1995.58) or indirectly from Thailand to Guangxi through Guangdong and then spread as a single lineage throughout the region. This result is in agreement with previous studies, suggesting that cross-border transmission of CRF01_AE between Vietnam and Guangxi occurs^[16]. In contrast, CRF08_BC strains, mainly originated in Yunnan, spread to Guangxi in 1992.29 (1989.05–1995.36) and formed the main epidemic strain among IDUs. CRF07_BC strains may originate in Yunnan and spread to Guangxi via Guangdong in 1995.23 (1993.64–1996.64). In contrast to previous studies^[18,29,30], our findings showed the source provinces of the majority of HIV-1 strains, except CRF08_BC, in recently infected IDUs were Guangdong, although we also observed some strains from Vietnam and Yunnan. Changes in the source provinces of these strains illustrate a shift in HIV-1 transmission drivers. They may have contributed to a new emergence of CRFs and URFs of HIV-1 and increased the complexity of the HIV epidemic among IDUs in Guangxi.

We identified three major CRFs (CRF08_BC, CRF01_AE, and CRF07_BC) in this study, which are consistent with previous results among IDUs in southwest provinces in China^[21] but were different from those in the Southeast coast and the Central Plains in China^[6,31]. Our study further confirmed that CRF08_BC is the dominant subtype present in IDUs in Yunnan and Guangxi and is rarely distributed in other places in the country. A previous study in China categorized CRF01_AE into seven independent lineages, namely, CRF01_AE Clusters 1–7^[32]. Our study identified four lineages (Clusters 1–4)^[23] among CRF01_AE strains circulating in Guangxi, but no strains from CRF01_AE Clusters 5–7 were found in this study. In addition, CRF01_AE Cluster 1 accounted for more than half of the total CRF01_AE strains, which is consistent with the findings of some recent studies^[33]. In our study, we observed a low prevalence of CRF07_BC strains and a high prevalence of CRF08_BC strains among IDUs in Guangxi, which is consistent with a systematic review reporting that CRF08_BC is the most prevalent genotype in the southwest region of China^[21]. Notably, CRF07_BC may be experiencing an increasing epidemic trend, especially among MSMs in Guangxi^[34], and we need to be alert to the possibility of a new wave of HIV infection in IDUs and MSMs driven by subtype CRF07_BC strains. Furthermore, we found 17 recombinant subtypes formed of inter-subtypes and inter-CRF, which was higher than those in previous studies^[6] and increased the complexity of HIV molecular patterns. This is also a major issue of a public health concern as novel strains may increase transmission, promote disease progression, or change resistance patterns to antiretroviral drugs^[35,36].

In the networks, we found that two large networks containing more than 100 IDUs infected with CRF01_AE or CRF_08BC were linked to heterosexuals from outside Guangxi and abroad, supporting the importance of identifying and controlling this “bridge population” to curb HIV transmission^[37]. The geographic dimension analysis (Figure 3) indicated that Guangdong and Yunnan were the major provinces for HIV-1 cross-regional transmission in Guangxi, accounting for 60.2% of all cases. High-speed railway and aviation may have facilitated the migration between these regions and promoted HIV-1 cross-regional transmission^[38]. Of course, there may be more fundamental social factors associated with the spatial mobility of HIV-1. As revealed by a previous study, the heterogeneity in the health status of urban dwellers could increase exposure and mobility and result in a high risk of disease transmission in large urban populations^[39]. Cities become incubators and provide conditions for outbreaks, and future migration and urbanization may strengthen disease transmission between cities^[39]. We also note that there was still an isolated cross-border transmission link to Vietnam. Although cross-border interventions have been effective in controlling the HIV epidemic among IDUs in recent years^[40], the risk of cross-border transmission still needs to be considered.

Although three common CRF strains in recently infected IDUs in Guangxi could be traced back to Vietnam, Guangdong, or Yunnan, the epidemic of HIV in Guangxi is mainly driven by local transmission. Specifically, in our study, more than 93% of HIV infections in the past 10 years were derived from HIV-1 strains already circulating in local Guangxi. Three common HIV CRF strains were most likely introduced to Guangxi IDUs in the early 1990s (1992–1995). Nonetheless, 6.4% IDUs were observed to be infected with HIV-1 from non-local circulating strains, which is comparable with the data for heterosexuals in Guangxi (5.9%)^[23]. Non-local circulating strains were mainly transmitted from Guangdong, consistent with our findings in geographic dimension analysis. Two putative reasons could explain the shift in HIV transmission sources among IDUs found in this study. First, in the mid-1990s, Guangxi experienced a high HIV prevalence among IDUs, which was derived from two emerging epidemics: one from Vietnam and another from Yunnan Province^[5]. In 2002, a Cross-Border Project for HIV prevention was well implemented among IDUs on both sides of the international border between Northern Vietnam and southern China^[40], playing an important role in reducing HIV epidemics for nearly a decade. Second, the high population mobility between Guangxi and Guangdong and high urbanized migration in Guangxi may contribute to the shift in HIV-1 transmission source. Accordingly, this study highlights the importance of understanding the role of external-province population mobility and intracity urbanized migration in HIV transmission.

An interesting finding of our study is that the primary HIV-1 transmission risk population among IDUs tended to be limited to PWID in recent years in Guangxi. Fortunately, these changes have not yet caused any rebound in HIV prevalence or incidence among IDUs. In fact, Guangxi suffered a massive HIV-1 epidemic among PWID between 1996 and 2003^[4]. Following the PWID epidemic in the 1990s, the predominant mode of HIV transmission shifted to heterosexuality through the bridge role of FSWs in HIV transmission from IDUs to the general population^[41]. In response to the HIV epidemic, the Guangxi government launched the Guangxi AIDS Conquering Project (GACP) in 2010^[42]. From 2009 to 2015, HIV/AIDS prevention efforts in Guangxi has been significantly strengthened, the scope of prevention services, MMT or needle exchange, peer education, and HIV testing has been expanded, and the awareness of HIV/AIDS knowledge among at-risk groups has been improved over time. These structural interventions have substantially altered HIV-related risk behavior and reduced HIV morbidity. Data from this study indicate that the local GACP initiative combined with the National AIDS Control Policy helps to stabilize HIV incidence in Guangxi.

Our study has several limitations. First, our study had limited samples and a convenience sampling method, both of which may lead to bias. Nevertheless, to the best of our knowledge, this study contains the largest number of sequences among IDUs in Guangxi to date to analyze the source of HIV strain transmission. Second, the missing variables of some classification indicators may have a potential impact on the analysis results. Third, most of the sequences from the other provinces were downloaded from the HIV database, as not all individuals living with HIV in other provinces have had their HIV-1 sequenced and submitted to the database. Last, we cannot identify the direct source of each transmission because there may be an intermediary between a case and its source. However, our analysis could identify closely linked infections.

In conclusion, we found a high diversity of HIV-1 strains among IDUs in Guangxi over the past decade and determined that HIV-1 transmission among IDUs in Guangxi mainly came from local strains that were already circulating. In addition to implementing traditional HIV prevention models, it is necessary to continue to strengthen the coverage of NSP and MMT. Furthermore, we have identified complicated sources of HIV transmission in Guangxi, suggesting appropriate measures to limit the spread of HIV along illicit drug-trafficking routes in the region.

The individual participant-level data that underlie the results reported in this article will be shared after de-identification (text, tables, figures, and appendices).