High Prevalence of HIV-1 Intersubtype B′/C Recombinants among Injecting Drug Users in Dehong, China

Objective To examine the distribution of HIV-1 genotypes among injecting drug users (IDUs) from Dehong, Yunnan province. Materials and Methods Blood samples from a total of 95 HIV-positive IDUs were retrospectively analyzed. Samples were collected between 2005 and 2009 from four cities in Dehong prefecture, western Yunnan province, the geographical origin of the HIV epidemic in China. HIV-1 gag, partial pol, vpr-env fragment, half-genome, or near-full-length sequences were analyzed to determine the HIV-1 genotypes of each subject. Results were compared with findings from past studies of IDUs in Dehong and in neighboring Myanmar. Results We observed a high prevalence of B′/C recombinants (82.4%) among IDUs in Dehong, the structural profiles of which do not match those previously reported in Dehong or in Myanmar. Furthermore, statistically significant differences in geographical and temporal distributions of HIV-1 genotypes were characterized by a predominance of HIV-1 B′/C recombinant forms among older subjects(p = 0.034), subjects from Longchuan district (p = 0.022), and subjects diagnosed between 2000 and 2004 (p = 0.004). Conclusions The increasing prevalence of multiple, new B′/C recombinant forms suggest that HIV-1 intersubtype recombination is substantial and ongoing in western Yunnan. This reflects the high-risk behavior of IDUs in this region and argues the need for stronger monitoring and prevention measures in Dehong and other high-prevalence areas around China.


Introduction
Yunnan province, located in southwestern China, was the first province to report evidence of an HIV epidemic and remains the most severely affected area in China today. The Dehong, Dai-Jingpo autonomous prefecture is one of eight minority autonomous prefectures in Yunnan province. According to China's 2010 national census, Dehong prefecture had a total population of approximately 1.2 million people with 49% belonging to minority groups, including Dai, Jingpo, Achang, Lisu, and Deang. Dehong encompasses an area of 11,526 square kilometers, divided into five districts: Luxi, Ruili, Longchuan, Yingjiang, and Lianghe ( Figure 1A). Except for Lianghe, each district shares a border with Myanmar. Travel, commerce, and intermarriage across the China-Myanmar border are common in this region.
Dehong has a long history of opiate trade, and the vast majority of illicit drugs in China are trafficked through Dehong from the ''golden triangle,'' which is comprised of large regions of Myanmar, Thailand, Laos, and Vietnam ( Figure 1A). In 1989, HIV-1 was first detected in injecting drug users (IDUs) in Ruili [1][2][3][4], and shortly thereafter in Longchuan and Luxi, all within Dehong prefecture [5][6][7]. These events marked the beginning of a true HIV epidemic in China. In 1992, a study examining HIV among IDUs in Dehong revealed that the prevalence of HIV infection had reached 81.8% in Ruili and 44.6% in Longchuan [6]. The spread of HIV in China was rapid and by 1995, outbreaks had been detected among IDUs in Sichuan, Xinjiang, and Guangxi provinces [8,9].
As the most recent molecular epidemiology study of HIV-1 in Dehong was performed before 2005 [11,15,16], the current status and trends of the HIV epidemic at the molecular level among IDUs in this region is largely unknown. Adjacent Myanmar is considered a geographical ''hot-spot'' of extensive HIV-1 subtype recombination with some 86.1% of HIV-positive IDUs in northern Myanmar carrying B9, C, and CRF01_AE recombinant strains of the virus. Although Dehong is expected to be in a similar situation because of the extensive trade and travel over the China-Myanmar boarder, the relationship between the HIV epidemics in Dehong and Myanmar remains unclear [20].
The primary aim of this study was to examine the distribution of HIV-1 genotypes among IDUs in Dehong, with secondary aims of investigating potential associations with demographic and epidemiological features of the study population and examining the current distribution in comparison to past findings.

Study population and blood sample collection
A total of 95 individuals were included in this study. Subjects were selected based on the following criteria: 1) self-reported acquisition of HIV through injecting drug use, 2) residence in one of the four Dehong districts that border Myanmar, and 3) enrollment between 2005 and 2009, at either the HIV/AIDS Care Center of Yunnan province or the Dehong Center for Disease Control and Prevention (Dehong CDC). Existing demographic information and clinical data as well as EDTA-3K-treated, anticoagulated whole blood samples were acquired for each subject. Samples were centrifuged at 300 g for 10 minutes and plasma was aliquoted and stored at 280uC. RNA extraction, amplification, and sequencing of HIV-1 gene fragments RNA was extracted from 280 ml of plasma using QIAampH Viral RNA Mini Kit (Qiagen, Germany) in a final elution volume of 60 ml. The sequences within the gag gene (HXB 2 790-2292), partial pol gene (HXB 2 2253-3318), and vpr-env region (HXB 2 5671-6449) fragments were amplified according to methods described in the literature [21][22][23]. In brief, HIV-1 gag, partial pol and vpr-env gene fragments were reverse-transcribed and amplified with SuperScript TM Polymerase One-Step RT-PCR System (Invitrogen, USA) and specific outer primers, respectively. Second round PCR amplification was performed using GoTaq DNA Polymerase (Promega, USA) and specific inner primers. PCR products were purified using QIAquick Gel Extraction Kit (Qiagen, Germany) and sequenced directly as described previously [21][22][23].
Single genome amplification and sequencing of HIV-1near-full-length genomes The 5-kb 59and39 half-genomes were amplified from RNA in plasma via RT-PCR with SuperScript TM III Reverse Transcriptase and Platinum Taq DNA Polymerase High Fidelity (Invitrogen, USA) as previously described [24]. Single genome amplification (SGA) and sequencing of HIV-1 DNA were performed in order to acquire the dominant single virus sequence from quasispecies [24]. Amplicons were sequenced directly by Beijing Genomics Institute (China) using internal walking primers.

Phylogenetic tree and recombination breakpoint analyses
Using the Los Alamos HIV Sequence Database (http://www. hiv.lanl.gov), all sequences were screened by the HIV BLAST tool to confirm lack of laboratory contamination and then aligned with HIV-1 reference strains. Alignment and manual editing were accomplished using Clustal X software (Version 2.0) and BioEdit software (Version 7.0; http://www.mbio.ncsu.edu/bioedit/ bioedit.html), respectively. Phylogenetic analyses were performed using the neighbor-joining method, based on the Kimura 2parameter distance matrix and a transition-to-transversion ratio of 2.0, using MEGA software (Version 5.0) [25]. Tree topology was tested by bootstrap analysis with 1,000 replicates. HIV-1 recombinant analysis was carried out using Simpolt (Version 3.5.1; http://sray.med.som.jhmi.edu/SCRoftware/simplot/) [26].The parameters of bootscan analyses were as follows: window size of 350 bp, step size of 50 bp for near-full-length and halfgenome sequences, window size of 200 bp, step size of 20 bp for gag, pol, and vpr-env gene fragments. The parameters of tree algorithm analysis were as follows: neighbor; distance model, Kimura; bootstrap replicate, 100; reference type, 50% consensus.

Statistical analysis
Categorical variables are presented as number and percent, while continuous variables are presented as mean and standard deviation (SD). Some continuous variables were divided into categories to examine threshold effects. Fisher's exact probability and one-way ANOVA analyses were used to compare demographic and clinical characteristics of the study population. HIV-1 genotype distributions across study sites, ethnic groups, and diagnosis times were also compared using Fisher's exact probability. P-values below 0.05 were considered significant. Missing data received no special treatment and were simply excluded from analyses. All statistical analyses were performed using SPSS software (Version 17.0).

Ethics statement
This study was reviewed and approved by the Medical Research Ethics Committee of No. 1 Hospital of China Medical University and written informed consent was obtained from all participants at the time blood samples were taken.

Nucleotide sequence accession numbers
The sequences reported in this article are available in GenBank under accession numbers of JX070462-JX070532 and JX070534-JX070556 for pol, KC189061-KC189103 for VPR to env, KC189104-KC189153 for gag and KC898975-KC899015 for near full length or half genome sequences.

Demographic and epidemiologic characteristics
The demographic and epidemiological characteristics of the 95 subjects included in the present study are shown in Table S1. The mean age of participants was 3567.

HIV-1 genotypes of IDUs in Dehong
From the 95 HIV-positive IDUs in this study, 33(34.7%) nearfull-length sequences, five (5.3%) 59 half-genomes, three(3.2%) 39 half-genomes, 88 (92.6%) gag gene sequences, 95 (100%) pol fragments, and 77 (81.1%) vpr-env fragments were successfully amplified and sequenced. Discrepancies in numbers of amplified fragments were due to the low volume of some specimens or low primer specificity. The genotype characterizations of HIV-1circulating among IDUs in Dehong were determined by phylogenetic and recombination analyses of the gene fragments of gag, pol, vpr-env, half-genome, and near-full-length sequences ( Figure S1 and Figure 2). The genotype characterizations of each subject as well as the genotypes of their specific gene fragments are listed in Table 1. Due to the complexity of the recombinants present in this study, the genotypes may have been different depending on the gene region analyzed. For example, the genotype B9/C accounted for 13.0% (10 of 77) of sequences in the vpr-env region, but 63.6% (56 of 88) of sequences in the gag gene region. Because of discrepancies such as these, only the 91 samples with near-full-length or half-genome sequences or at least 2 sequences of pol, gag and vpr-env fragments were used in the following analyses on genotype-related factors.

Factors associated with HIV-1 genotypes
Subjects were grouped into one of four categories (B9, B9/C, C, and Other) to represent their HIV-1 genotypes ( Table 1).The genotype distribution of the 91 IDUs used for genotype-related factors analyses were as follows: 75 (82.4%) B9/C recombinants, 5 (5.5%) B9 genotypes, five (5.5%) C genotypes, and six (6.6%) other genotypes, including three that were B9/C/CRF01_AE, two that were C/CRF01_AE, and one that was pure CRF01_AE. The relationships of HIV-1 genotypes to the demographic and epidemiologic features of study population are summarized in Table 2. A statistically significant association was observed between HIV-1 genotype and age (p = 0.034), district (p = 0.022), and year of diagnosis (p = 0.004), suggesting both geographical and temporal differences in the distributions of the different HIV-1 genotypes. No other statistically significant associations were detected.

Geographical distribution of HIV-1 genotypes
As shown in Table 2 and Figure 1B, the prevalence of B9/C recombinants was more than 70% in Longchuan, Luxi, and Yingjiang, with the highest prevalence observed in Longchuan (45 of 48, 93.8%),where the prevalence of subtype C was found to be the lowest (0 of 48) among the four districts. By contrast, the lowest prevalence of B9/C recombinants was found in Ruili (15 of 23, 65.2%), where pure subtype B9 and subtype C was at its highest prevalence (6 of 23, 26.0%). This higher prevalence of pure subtype B9 and subtype C in Ruili relative to other districts was statistically significant (p = 0.022). It is important to note that 5 of the 6 pure subtype B9 or subtype C genotype HIV-positive subjects were diagnosed before 1995 and were between 40 and 48 years old.
In order to examine potential geographical linkage with nearby regions, we compared the HIV-1 genotypes of IDUs in Dehong with data from previous HIV-1 genotyping studies on IDUs in Myanmar. As shown in Table 3, significantly more B9/C recombinants and fewer other genotypes were found in the Dehong study population when compared with findings from both northern and central Myanmar (p,0.001) [20,27]. In order to elucidate the phylogenetic relationships between HIV-1 strains among IDUs in Dehong and Myanmar, neighbor-joining tree analysis was conducted based on near-full-length sequences and gene fragments of gag, pol, vpr-env regions from this study as well as two previous studies (Figure 3 and Figure S2). Although various HIV-1 recombinant forms, especially B9/C recombinants, were detected in both Dehong and Myanmar, none of the recombinant sequences from Dehong and Myanmar clustered together with a bootstrap value above 70 in near-full-length neighbor-joining tree analysis except for a cluster that included three sequences from Dehong IDUs, three sequences from Myanmar IDUs, and one sequence from a Myanmar patient who acquired HIV via sexual contact (bootstrap value 100; Figure 3A). However, breakpoint analysis demonstrated these sequences had different recombination patterns or no breakpoints ( Figure 3B).Two other clusters of Dehong and Myanmar sequences were observed in pol neighborjoining tree analysis (bootstrap values of 92 and 98). However,   vpr-env (5671-6449) these sequences did not cluster together in the gag and vpr-env neighbor-joining tree analyses ( Figure S2), suggesting that although the recombinants from Dehong and Myanmar are closely related, the recombination events that gave rise to these strains occurred independently in Dehong and Myanmar.

Temporal distribution of HIV-1genotypes
In our samples, we observed a trend of increasing prevalence of HIV-1 B9/C recombinants in IDUs from Dehong who were diagnosed at different periods ( Table 2). In those diagnosed before 1995, 73.3% carried B9/C recombinant genotypes, whereas in those diagnosed between 1996 and 1999, the proportion had increased to 88.9%, and then between 2000 and 2004, had increased further to 94.3%. However, after 2005, the prevalence of B9/C recombinants dropped dramatically to 44.4%. Over this same period, from prior to 1995 to after 2005, a few pure subtype B9 and pure subtype C genotypes were detected, though these were very rare.
In order to more closely examine temporal changes in genotype distribution, we compared our results with data from previous HIV-1 genotyping studies among Dehong IDUs. As summarized in Table 3, HIV-1 subtype B9 was the predominant form (49 of 54, 90.7%) in samples collected from 1996 to 1998 [15], while B9/C recombinants were less common (4 of 54, 7.4%). In contrast, studies performed on samples collected from 2000 to 2001 [11] and in 2003 [16] document the prevalence of B9/C recombinants as being greater than 60% and the prevalence of subtype B9 at less than 30%. Taken together, our results and those of previous studies suggest that while B9/C recombinant forms are arising continuously and have become the dominant HIV-1 subtype in Dehong prefecture, pure subtype B9 and pure subtype C genotypes have not yet disappeared from this population.

Evolution of HIV-1 recombinants in Dehong
To further explore the evolutionary relationships of B9/C recombinants in Dehong, and to investigate the possibility of second generation recombination, we compared the recombination patterns of 33 near-full-length sequences from IDUs diagnosed at different time periods. In the phylogenetic tree generated using near-full-length sequences, 25sequences were clearly grouped into nine clusters (bootstrap value .80; Figure 2A). In addition, seven B9/C strains and one pure B9 strain was observed. The sizes of the nine clusters ranged from two to five sequences. The results of our recombination analysis reinforced the findings of our phylogenetic analysis: at least 16 different recombination patterns were observed, none of which match those previous reported in Dehong ( Figure 2B and Figure 2C). In general, while the 25 sequences in clusters have recombinants from B9 and C subtypes, it is difficult to resolve the evolutionary relationships between them. However, it is interesting to note that the 8 sequences in clusters 4, 5, and 9 all had a B9 fragment in the env-nef region and all shared the same breakpoints ( Figure 2B and Figure 4A). The sub-region tree suggested that sequences in clusters 4, 5, and 9 had a common subtype C backbone and B9 In the recombination analyses, a subtype C strain from India (95IN21068) (pink) and a subtype B9 strain from Yunnan (RL42) (blue) were used as subtype references, and a CRF01_AE strain from Thailand (CM240) (green) was used as an out-group control. doi:10.1371/journal.pone.0065337.g003 fragment insertion in the env-nef region ( Figure 4B). However, the sequences in clusters 4 and 9 contained more B9 fragment insertions in the 59 half-genome, compared with those in cluster 5, implying the possibility of second generation recombination. In addition, the B9 fragments of clusters 4 and 9 in 59 half-genome clustered together, implying a homologous relationship  ( Figure 4B).The recombination maps of the 26 near-full-length sequences and 6 half-genomes from patients that have the diagnoses time are presented in Figure 5. These data show that more pure B9 and pure C sequences are observed in earlier groups, while more B9/C recombinant forms are observed in later groups. However, in the four sequences from subjects diagnosed between 2008 and 2009, two pure subtype C, one B9/C recombinant, and one CRF_01_AE/B9/C recombinant were observed.

Discussion
In the present retrospective study, using blood samples collected from 95 HIV-positive IDUs from Dehong between 2005 and 2009, we determined HIV-1 gag gene, partial pol gene, vpr-env region, half-genome, and near-full-length nucleotide sequences and analyzed their phylogenetic relationships and recombinant structures. These were compared to the demographic and epidemiological data of our study population as well as to results from prior HIV-1 molecular epidemiology studies both in Dehong and in Myanmar.
The main finding of our study was the significant predominance of HIV-1 B9/C recombinant forms among Dehong IDUs (Table 1 and Figure 1B), which have structures that have not previously been documented in the literature (Figure 2). The different B9/C recombinant forms we have found may have been generated via recombination events occurring not only between subtype B9 and subtype C forms, but also between different B9/C forms via second generation recombination events. The phylogenetic and recombination analysis on the near-full-length sequences supported the idea that sequences in clusters 4 and 9 may be second generation recombinant forms that have arisen from B9 and B9/C recombinants in cluster 5 ( Figure 2). We believe this may be further supported when the year of diagnosis is considered. Three of the four subjects each in clusters 4 and 9 were diagnosed in 1996 (the fourth was diagnosed in 2004). Unfortunately, diagnosis year data was unavailable for the four sequences in cluster 5, so we cannot draw conclusions about these potential second generation recombinants based on subject diagnosis year. However, our phylogenetic analyses showed likely evolutionary relationships among some B9/C recombinants. Furthermore, the increasing complexity Figure 5. The bootscanning plots of 26 near-full-length sequences and 6 half-genome sequences of HIV-1 strains isolated from Dehong IDUs with diagnosis year in chronological order. In the recombination analyses, a subtype C strain from India (95IN21068) (pink) and a subtype B9 strain from Yunnan (RL42) (blue) was used as subtype references, and a CRF01_AE strain from Thailand (CM240) (green) was used as an out-group control. doi:10.1371/journal.pone.0065337.g005 of recombination patterns suggests that there is a high frequency of HIV-1 recombination events within this population and new recombinant forms are arising continuously ( Figure 5). This phenomenon has been documented in neighboring Myanmar, where 86.1% of samples in one study were found to be recombinants [27]. The co-circulation of multiple lineages of HIV-1 strains in pre-existing transmission networks with high-risk behaviors (like that of IDUs) could easily fuel the generation of a variety of new recombinant strains in China.
This idea is further supported by evaluating our findings in the context of previous molecular epidemiological studies of HIV-1 recombinant forms in Dehong prefecture. The high prevalence of B9/C recombinants in our sample is consistent with previous observations by Yang et al. for samples collected between 2000 and 2001 [11], and Zhang et al. for samples collected in 2003 [16], but is different from samples collected by Qiu et al. between 1996 and1998 [11,15,16], in which HIV-1 B9 subtype virus was most common. While this makes sense in light of the fact that B9/C recombinant forms first became common in Dehong in 1996, shortly after the pure HIV-1 subtype C virus was introduced [11][12][13][14][15][16]18]. These data taken together, point to an increase in the prevalence of B9/C recombinants over time (Table 3).
Surprisingly however, two of the four cases diagnosed between 2008 and 2009 (with seronegative records within six months), were identified as having pure subtype C strains, seemingly contrary to the trend of increasing B9/C recombinant prevalence. It is important to note that in recent years, many HIV prevention policies have been implemented by the Chinese government. In Dehong, needle exchange programs began in 2003, and the first methadone maintenance treatment clinics were established in 2005. These newly-infected subjects carrying pure subtype C viral strains are therefore not infected with multiple strains of HIV-1.This may be an indication that the HIV prevention policies have achieved some success, even in areas with high rates of recombination such as Dehong.
Our results regarding the geographical distribution of HIV-1 genotypes are notable and in line with our analysis of the HIV-1 genotype temporal distribution in Dehong (Table 2and Figure 5). While the prevalence of HIV-1 B9/C recombinants is relatively high in all four districts, significantly higher prevalence of pure subtypes B9 and C were observed in Ruili, relative to other regions. A possible explanation for this finding is that subjects from Ruili were infected with HIV in the early-and mid-1990s when fewer B9/C strains were present in the emerging epidemic in Dehong (5 of 23, 21.7% of Ruili subjects were diagnosed before 1995, data not shown).
Another significant finding is the relationship of the HIV epidemics in Dehong and neighboring Myanmar. Dehong is the most active channel for the trafficking of illicit drugs from Myanmar and other golden triangle countries into China. A high prevalence of diverse forms of HIV-1 B9, C and CRF01_AEintersubtyperecombinantshas been reported in Myanmar [27,28]. With the most recent report from northern Myanmar finding that 86.1% of HIV-positive IDUs carry intersubtype recombinant forms of the virus [20]. The recombinant forms identified in this report include all four possible strains formed from recombination events between pure subtype B9, pure subtype C, and CRF01_AE. In addition, most recombinants have distinguishable chimeric patterns, forming 64 URFs that are dissimilar to previously identified CRFs and other URFs in Asia. In the present study, we detected a similar prevalence of intersubtype recombinant forms in Dehong (80 of 91, 87.9%), which prompted us to consider the possibility of a close relationship between the HIV epidemics in the two regions. However, comparison of the phylogenetic relationship among recombinants in Dehong and Myanmar revealed that none of the available sequences shared the same recombination patterns. Taken together, these data suggest that Dehong and Myanmar are both experiencing high frequency and complex HIV-1 intersubtype recombination events, but that they are largely independently of each other.
Our study had several limitations. Firstly, comparison of results from this study to previous molecular epidemiological observations is difficult due to inter-study variability of the gene regions analyzed. For example, Qiu et al. used the C2V3 region of env gene to determine genotype [15]. Because breakpoints are more commonly documented in the pol gene [8][9][10][11], it could be that the prevalence rates of HIV-1 B9/C recombinant genotypes in the Qiu et al. study were underestimated. Secondly, because subjects who were diagnosed in earlier years could have had a greater probability of being infected with HIV multiple times (compared to subjects diagnosed in later years), their samples may have presented with more recombinants. Thirdly, our dataset, while much larger in sample size than previous studies, did have some missing values. This could have resulted in the under-or overestimation of the association of HIV-1 B9/C recombinant genotypes with demographic or epidemiologic factors.
In conclusion, our study documents a high and increasing prevalence of HIV-1 B9/C recombinants among Dehong IDUs since 1996, shortly after the subtype C virus was introduced in the region. Further, the increasing predominance of multiple, new, forms of unique B9/C recombinant strains suggests that HIV-1 intersubtype recombination is extensive and ongoing among Dehong IDUs. These findings highlight the urgency of strengthening monitoring efforts and implementing effective measures to reduce transmission in this population. It also underscores the need for further research on the rate of spread of these forms outside of Dehong and the potentially altered biology of these new recombinant forms. The increasing mobility of people across local and international borders, presents a high potential for the rapid evolution of an HIV strain with exaggerated transmission efficiency, virulence, or drug resistance, which is a global problem of critical importance. Figure S1 Neighbor-joining trees of gag, pol and vpr to env fragments from Dehong IDUs. Neighbor-joining trees were built based on the sequences of 88 gag gene (A), 95 pol gene (B) and 77 vpr to env gene fragments. The subtype references sequences from the Los Alamos HIV Sequence Database were included as the reference sequences. The stability of the nodes was assessed by bootstrap analyses with 1000 replications. Only bootstrap values of more than 70 are shown at the corresponding nodes. The sequences from IDUs in Dehong are labeled by red solid circles, and previous published sequences from Dehong are labeled by black triangles. The clusters that including sequences from Dehong both in this study and previous published sequences were showed in red line.