Geographic Transmission and Epidemic History of HIV-1 CRF01_AE, CRF07_BC, and HCV Subtype-6w among Taiwanese Persons Who Inject Drugs

Persons who inject drugs (PWID) and their risk-related behaviors (e.g., unprotected sex and sharing needles/syringes/other injection equipment) have caused severe public health problems, especially in the rapid spread of HIV-1 and HCV. Here, we reconstructed the epidemic history of HIV-1 circulating recombinant form (CRF) 01_AE, CRF07_BC, and HCV subtype-6w among Taiwanese PWID. The timescales were estimated using phylogenetic and Bayesian coalescent analyses. The results revealed that CRF01_AE started to circulate in the Taiwanese PWID population in central Taiwan at 1992.5 (95% credible region: 1988.8–1995.9) and spread to other regions of Taiwan, while CRF07_BC was first identified in southern Taiwan at 2000.0 (95% CR: 1997.8–2002.2) and then spread northward to central-northern Taiwan. All HCV-6 strains were from Asia (that is, China, Myanmar, Taiwan, and Vietnam) and originated in 1928.1 (95% CR: 1890.2–1966.0). Furthermore, subtype-6w isolates from different regions of Taiwan appeared to share a common source that existed in the mid-1990s (95% CR: 1985.0–2001.8) or thereabouts. The routes of drug trafficking and the resulting high prevalence of HIV-1/HCV co-infections among PWID might have contributed to the virus transmission and promoted its spread worldwide. Long-term monitoring and policy implementation in at-risk populations would be useful for disease control.


Introduction
In 2021, 38.4 million people globally were estimated to have lived with HIV and the worst region is still sub-Saharan Africa [1]. A systematic review revealed that of the 16 million people who inject drugs (PWID), one-fourth lived in East and Southeast Asia [2]. PWID and their risk-related behaviors (e.g., unprotected sex and sharing syringes/heroin solutions/other injection equipment) have caused server public health problems, especially regarding the rapid spread of HIV-1 and HCV. According to the distribution data for the HIV-1 subtypes and the circulating recombinant forms (CRFs) in Asia, the HIV-1 subtype B' and CRF01_AE are the main subtypes among PWID in Thailand [3][4][5], Myanmar [6], and Vietnam [7], while CRF07_BC is the predominant strain among PWID in China and Taiwan [8][9][10][11][12]. The strain CRF07_BC is derived from the Thai subtype B' and Indian subtype C lineages [13]. Since 1997, CRF07_BC has been isolated from PWID in several

Subjects
The research procedures for the current study are shown in Figure A1 analysis. All Taiwanese PWID were obtained by direct sequencing in our laboratory and those retrieved from GenBank were listed in Table A2.
Sociodemographic data and information on the types of illegal drugs used, history of drug abuse, risk factors associated with HIV-1 transmission, and years of the first HIV/HCV positive diagnosis were collected using a self-administered questionnaire. Peripheral blood samples were collected to allow analysis of virus genotype. Informed consent was obtained from all participants. Our research protocol was approved by the prisons and detentions administration system, as well as the Institutional Review Board of the National Yang-Ming University, Taiwan.

HIV-1 and HCV Subtyping
Viral RNA was extracted from plasma samples using the QIAamp Viral RNA mini kit (QIAGEN, Hilden, Germany). Random primer (Promega) was used in reverse transcription to generate cDNA for reverse transcriptase-polymerase chain reaction (RT-PCR). Anti-HCV antibodies from serum samples were detected using an enzyme immunoassay system (Murex 3nd, Abbott Laboratories, North Chicago, IL, USA). Specimens determined with anti-HCV antibodies or confirmed as HIV-1 positive were further analyzed. The genotypes/subtypes of HCV and HIV-1 infections were determined according to the methods described previously [10,25,32]. A set of primers, OF9-2 (forward) 5 -CGACATTACGCAGAAGTTGCCC-3 and OR9 (reverse) 5 -AGTGTTGCTTAAGGCCTCCTGC-3 , were used to amplification of the HCV NS5B gene near full-length. Proviral nucleotide sequences were obtained by direct sequencing of PCR products using a DNA analyzer (ABI 3730, Applied Biosystems, Foster City, CA, USA).

Phylogenetic Analysis
Sequence alignment analysis with various reference strains from the Los Alamos HIV-1 database (https://www.hiv.lanl.gov/content/index, accessed on 28 May 2021) and the HCV database (https://hcv.lanl.gov/content/index, accessed on 8 July 2021) was performed using the BioEdit v7.2.6.1 program [33]. The MEGA X program [34] was used to find the best fit nucleotide substitution model and to construct phylogenetic trees using neighbor-joining (NJ) and maximum likelihood (ML) methods. For example, taking the HIV-1 env gene, the substitution model GTR + G was incorporated into the ML method, while TN93 + G (GTR and HKY models are not available here) was used to calculate the evolutionary distance for the NJ tree followed by bootstrap analysis with 1000 replicates [35]. Considering the best-fit models for the HCV NS5B gene, the substitution models for both the ML and NJ tree were K2 + G + I. At least four strains of all subtypes were used as reference sequences and isolates from Asian PWID were included for phylogenetic analysis (see Figure A2). Bootstrap values (≥70%) were used as an indicator of the significance of the clusters.

Nucleotide Sequence Accession Numbers
The HIV-1 env sequences (OM287868-OM287928) and the HCV NS5B sequences (OM287929-OM287937) were obtained from the current study and deposited in GenBank.

Bayesian Coalescent Inference
Evolutionary rates were obtained using the Bayesian Markov chain Monte Carlo (MCMC) approach implemented in BEAST v2.5.1 [36]. General time-reversible (GTR) [37][38][39] substitution models with gamma-distributed among-site rate variation involving six categories [40] were used to estimate evolutionary rates and construct tree topologies. Constantly sized, exponentially growing, and Bayesian skyline coalescent models were used for each case [41] and each MCMC chain was run for at least 10,000,000 states and sampled in every 1000 states. Posterior probability densities were calculated, and the convergence of the chains was verified using the Tracer v1.7.1 [42] with 10% of each chain discarded as burn-in.

Statistical Analysis
The Pearson χ 2 test and Fisher's exact test were performed in univariate analysis of demographic data. The difference between groups with a p-value < 0.05 was considered statistically significant. The p-values were two-tailed and unadjusted for multiple comparisons.

Results
3.1. Geographical Distribution of HIV-1 CRF01_AE, CRF07_BC, and HCV Subtype-6w among Taiwanese PWID From 2004 to 2005, almost all Taiwanese HIV-1 positive PWID were infected with CRF07_BC. However, we found another small-scale outbreak strain that circulated in this population in central-north Taiwan. These were judged according to their time of crime and the place of sentences.
As shown in Figure A3, the distribution of Taiwanese PWID infected with HIV-1 CRF01_AE (n = 24), HIV-1 CRF07_BC (n = 982), and HCV subtype 6w (n = 9) during 2004-2008. The dates of HIV diagnosis and sample collection for most CRF01_AE infections were mainly in 2005 ( Figure A2 and Table A2). The CRF07_BC sequences of Taiwanese PWID were grouped into several distinct phylogenetic clusters based on collection places [11] (details shown in Table A2). Based on the dates of HIV diagnosis, our data implied that CRF01_AE started to circulate in the Taiwanese PWID population in Central Taiwan and then spread to other regions of the island. In contrast, CRF07_BC first appeared in the south and moved northward to expand to central-north Taiwan ( Figure A2 and Table 1).
The nine Taiwanese PWID infected with HCV-6w were identified when serving their prison sentences. Among them, five were HCV mono-infections and the other four cases were HIV/HCV coinfections. Six of the cases were from northern Taiwan (Taipei Detention Center, Sindian Drug Abuse Treatment Center and Taoyuan Woman's Prison), two were from central Taiwan (Taichung Prison and Yunlin Second Prison), and one was from southern Taiwan (Kaohsiung Prison) ( Figure A3).

An Estimated Timescale of the Spread of CRF01_AE and CRF07_BC among Taiwanese PWID
When estimating the time scale of the spread of HIV-1 CRF01_AE and CRF07_BC in Asia, we adopted the data set based on GTR + Γ 6 constant model to pinpoint the time of the most recent common ancestor (tMRCA) of the HIV-1 strains circulating in this area. Compared to the results of the three models, similar conclusions could be reached (Table 1). After systematic analyses, we followed the likelihood of constant size, exponential growth, and the Bayesian skyline model (CRF01_AE: −15,559.7307, −15,548.8808, and −15,719.0827; CRF07_BC: −4939.2117, −4933.7994, and −4957.3370) and found that the exponential growth model was the best to present its transmission.
The estimated phylogeny using the env gene showed that all Taiwanese CRF01_AE and reference strains formed a single clade. However, there are two different risk groups in Taiwan [11], namely PWID and other sexual groups (e.g., homo, hetero, and bisexuals). The former population that contained sequences from Taipei Figure 1a). A comparison with our previous findings [11] suggested that CRF01_AE was introduced into Taiwanese PWID through unprotected sex and then caused a local epidemic among PWID through the exchange of injection equipment.  Figure 1a). A comparison with our previous findings [11] suggested that CRF01_AE was introduced into Taiwanese PWID through unprotected sex and then caused a local epidemic among PWID through the exchange of injection equipment. 2) in Northern, Central, and Southern Taiwan, respectively. The CRF07_BC strains from different regions of Taiwan seem to share a common source that existed in 2000 or thereabouts (95% CR: 1997. 8-2001.9) and was part of the Southern Taiwan PWID. This suggests that southern Taiwan was the entry site for CRF07_ BC (Table  1 and Figure 1b).

An Estimated Timescale of the Spread of HCV Subtype-6w among Taiwanese PWID
Similarly to the estimation of the timescale of HIV-1 spread, we adopted the data set based on GTR+Γ6+I constant model to pinpoint the tMRCAs of HCV-6 and found that the exponential growth model (likelihood in CS: −9744.679, EG: −9717.44, and BS: −9868.01) is the best way to present its transmission.
Phylogeny analysis using the NS5B gene showed that all Asian PWID and reference strains formed a single clade. As summarized in Table 2, all genotype-6 strains from Asia (i.e., China, Myanmar, Taiwan, and Vietnam) were rooted in 1928.1 (95% CR: 1890. .0). Subtypes 6a, -6n, and -6w had existed in the Taiwanese PWID population (Table  A1). Taking subtype-6a for example, it was introduced into Vietnam at 1993.5 (95% CR: 1977.5-2001.3) and later into China at 1994.5 (95% CR: 1988.9-2000.9). Subtype-6n was initially introduced to China (1987China ( .8, 1952.0-2005.0) and then spread to Myanmar (1990Myanmar ( .4, 1954Myanmar ( .7-2007. It is noteworthy that this strain was found to originate in Yunnan (1987. 8,1953.1-2007.0) and spread eastward to Suzhou, Zhenjiang, and Jiangsu in the early and mid-2000s. Furthermore, subtype-6w isolates from different regions in Taiwan seem to share a common source that existed in mid-1990 (95% CR: 1985.0-2001.8) or thereabouts ( Table 2 and Figure 1c). 2) in Northern, Central, and Southern Taiwan, respectively. The CRF07_BC strains from different regions of Taiwan seem to share a common source that existed in 2000 or thereabouts (95% CR: 1997. 8-2001.9) and was part of the Southern Taiwan PWID. This suggests that southern Taiwan was the entry site for CRF07_ BC (Table 1 and Figure 1b).

An Estimated Timescale of the Spread of HCV Subtype-6w among Taiwanese PWID
Similarly to the estimation of the timescale of HIV-1 spread, we adopted the data set based on GTR + Γ 6 + I constant model to pinpoint the tMRCAs of HCV-6 and found that the exponential growth model (likelihood in CS: −9744.679, EG: −9717.44, and BS: −9868.01) is the best way to present its transmission.

Discussion
Takebe et al. reported that CRF07_BC strains from different regions in China (including Xinjiang, Liaoning, and probably Guangdong and Sichuan) were likely to share a common ancestor that existed in Yunnan province around 1993 (95% CR: 1991.2-1995.2; gag) [13,31]. This suggests that CRF07_BC spreads almost simultaneously to various regions of China [13,31]. Furthermore, CRF07_BC also spread to Taiwan from the South around 1999.7 (95% CR: 1998.  [10,13,30], resulting in a major HIV epidemic among PWID in Taiwan [10][11][12][13]31]. The dissemination routes of CRF07_BC in China and Taiwan were those reported in previous studies [11,13,30,31]. To compare the main differences between the use of the date of sample collection versus the date of HIV-1 diagnosis to estimate the time of the emergence of the CRF01_AE and CRF07_BC strains, and to consolidate the integrity of our data, we added more sequences from West Taiwan (e.g., Sindian, Taoyuan, Taichung, Yunlin, and Kaohsiung) in the analysis. As we all know that CRF07_BC circulated in southern Taiwan first, even adding the sequences from the most south area (i.e., Kaohsiung), the tMRCAs of CF07_BC among Taiwanese PWID were behind the estimates as previously reported [13,30,31]. The data obtained using the date of sample collection are more accurate than those using the date of diagnosis. This finding showed why it is necessary to use the correct date to estimate the time of emergence of HIV-1 subtypes or CRFs. Furthermore, our results revealed that the estimated introduction time of CRF01_AE in Taiwan PWID (1992 later) was earlier than that of CRF07_BC (1999.9), and because of the less aggressiveness of CRF01_AE [43], it only caused a local epidemic initially.
The estimated prevalence of HCV-6 in some regions of Southeast Asia, especially among patients with PWID and major thalassemia, is as high as 50% [44]. Injecting drug abuse is possibly responsible for the high frequency of this genotype in certain parts of Asia. HCV-6 has considerable genetic diversity with 23 subtypes (a-w). HCV-6 infected with HCV-6 respond better to interferon-based therapy compared to genotype 1, although the clinical characteristics and side effect profiles in patients are similar between HCV-6 and other genotypes [44]. Our study showed that HCV-6 was as common as genotype 1 (34.7% vs. 43.5%, Table S1) in the Taiwanese PWID population. According to a largescale survey on the seroprevalence of HCV in Taiwan [45], the prevalence of injecting drug abuse and incomplete disinfection of medical utensils would cause a small-scale outbreak of HCV in local areas. Furthermore, residents have a higher prevalence of HCV when they were born in an earlier cohort [45]. As shown in Table A1, a crosssectional study with 624 PWID recruited in Taiwan was conducted in 2007-2008. The overall prevalence of HIV and HCV infection was 44.1% (275/624) and 80.4% (502/624), respectively. The prevalence of HCV mono-infection and HIV/HCV co-infection was 36.4% (227/624) and 44.1% (275/624), respectively. The issues of HCV prevention include the following: to prevent healthy people from being contaminated with infected blood and to avoid reinfection with HCV in cured cases. For those who have been cured and non-infected, regular screening tests are encouraged. Through a series of analyses, our findings appear to support the hypothesis that HCV-6 originated in Southeast Asia (Table 2). HCV-6 is highly divergent from other genotypes and has distinct genetic differences from other strains, suggesting that there may be unclassified subtypes in Asia. Therefore, the accumulation of such genetic heterogeneity suggests that this genotype has circulated, adapted, and evolved in this area for a long period.
There are several limitations to this study. First, all Asian isolates (e.g., China, Malaysia, Myanmar, and Vietnam) were restricted and obtained from the NCBI website. Second, some Asian isolates were excluded from the evolutionary analysis because the sequences were too short or contained missing sequences. Despite the limitations, this study sheds light on the routes of drug trafficking and the resulting high prevalence of HIV-1/HCV coinfections among PWID that could have contributed to regional and global transmissions. In conclusion, for the first time, we report 'the time of emergence of common HCV and HIV-1 strains among Taiwanese PWID' and provide a comprehensive profile suggesting the initial circulation of CRF07_BC in southern Taiwan before spreading to other regions of Taiwan. Furthermore, the importance of using the date of sample collection versus the date of HIV-1 diagnosis was also highlighted when estimating the time of the emergence of the CRF01_AE and CRF07_BC strains. Long-term monitoring and implementation in the population at risk would be useful for disease control. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Also written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
The data presented in this study are available in [Section 2.4 and Table A2]. Part of this paper was presented at the 8th Second Member Conference and Academic Symposium held in Taiwan in 2020.

Acknowledgments:
We thank the PWID who participated in this study, the peer educators and social workers from the AIDS Prevention and Research Center who helped collect the questionnaires, and the staff of the Genomic Research Center of the National Yang Ming Chiao Tung University, as well as the health clinics of the various prisons and detention centers for their administrative support and technical assistance. The Sequencing Core Facility is supported by the National Research Program for Genomic Medicine (NRPGM) of the National Science Council. Most importantly, we would like to thank Wing-Wai Wong, who provided clinical samples but unfortunately passed away in October 2018, for his contribution to this study.

Conflicts of Interest:
All authors of this article do not have commercial or other associations that might pose a conflict of interest. In addition, the funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.              Note: a The lengths of the env fragment in CRF01_AE and CRF07_BC were up to 537 and 552 bp, respectively. The length of the query shorter than the standard was labeled with a different color, such as 'light gray, <500 bp' and 'dark gray, <300 bp'.