In this study, 3,596 individuals with newly confirmed HIV infections from 2020 to 2023 and 164 individuals with ART failure from 2021 to 2023 were enrolled from Hangzhou, Zhejiang Province, China. Blood samples and socio-demographic characteristics, including sex, age, ethnicity, education, current place of residence, marital status, infection route, and high-risk sexual behavior, were collected. All participants provided informed consent, and the study was approved by the Medical Ethics Committee of the Hangzhou Municipal Center for Disease Control and Prevention.

Plasma was separated from whole blood at 3,000 × g for 10 min and then stored at −80°C. HIV-1 RNA was extracted from the plasma using the QIAamp Viral RNA Mini Kit (QIAGEN, Germany). The pol gene fragments (1,316 bp, HXB2: 2147–3,462) were amplified using reverse transcription polymerase chain reaction (RT–PCR) and nested polymerase chain reaction (PCR). The amplification products from the PCR-positive samples were sent to Hangzhou TsingKe Bio-Tech Co. for purification and sequencing. The obtained sequence fragments were assembled using Sequencher 5.4.6 (Gene Codes, Ann Arbor, MI, United States).

The sequences obtained were aligned with the HIV reference sequences downloaded from the Los Alamos National Laboratory HIV Sequence Database and corrected using MEGA v11.0.13. Subtypes were identified using HIV BLAST.¹ The sequences were then uploaded to the Stanford HIV Drug Resistance Database² for the analysis of resistance mutation sites.

The HIV-1 molecular transmission network was established using the HIV-TRACE algorithm (Kosakovsky Pond et al., 2018). Pairwise genetic distances of the sequences were calculated with the TN93 model using HyPhy (Pond et al., 2005). The HIV molecular transmission network was constructed with a 0.5% distance threshold, which is recommended by the CDC in the United States and China (CDC, 2018; Oster et al., 2018). The network was visualized in Cytoscape v3.10.1.

Forty-four CRF07_BC reference sequences from other regions in China were downloaded from GenBank based on their high homology (>98.0%) with the participants in drug-resistant transmission clusters. BEAST v.1.10.4 under an uncorrelated relaxed clock model, the GTR + G + I nucleotide substitution model, and a Bayesian skyline plot demographic model were used to perform Bayesian evolutionary analysis (Suchard et al., 2018). The BEAST analysis was performed using Markov Chain Monte Carlo (MCMC) runs of 100 million generations, with samples taken every 10,000 steps. The Bayesian MCMC output was analyzed using Tracer v1.7.2 (Rambaut et al., 2018). Maximum clade credibility (MCC) trees were generated using TreeAnnotator v1.10.4 and visually edited in FigTree v1.4.4.

A chi-squared test, Fisher’s exact test, and two-way ANOVA were performed in GraphPad Prism 9. Turkey multiple comparisons test was performed after the two-way ANOVA. The Chi-square test was employed as univariable regression. Multivariable logistic regression was performed in SPSS 25. The odds ratios (OR), adjusted odds ratios (aOR), and 95% confidence intervals (CI) were calculated.

Out of 4,595 plasma samples from newly confirmed HIV infections and 211 plasma samples from individuals with ART failure, 78.2% (3,596/4595) and 77.8% (164/211) were successfully sequenced, respectively. The participant counts of newly confirmed HIV infections and individuals with ART failure in the following analyzed datasets were N1 = 3,596 and N2 = 164, respectively (Table 1). There was no significant difference between the participants and non-participants in our analysis.

Among the newly confirmed HIV infections, the majority were male (90.5%, 3,254/3,596), younger individuals (78.0%, 2,804/3,596), and of Han ethnicity (95.3%, 3,426/3,596). Over half had at least a high school education (61.7%, 2,217/3,596), were unmarried (57.1%, 2,053/3,596), acquired HIV through homosexual contact (65.7%, 2,327/3,596), and had a CD4+ T cell count between 200 and 500 cells/μL (53.2%, 1,913/3,596). The most common HIV subtype was CRF07_BC (45.8%, 1,648/3,596), followed by CRF01_AE (33.9%, 1,219/3,596), CRF08_BC (5.2%, 188/3,596), CRF55_01B (4.4%, 159/3,596), URF (CRF07_BC/CRF01_AE; 3.1%, 111/3,596), and the subtype B (2.3%, 82/3,596). The individuals with ART failure exhibited similar proportions in terms of sex, age, ethnicity, education level, marital status, and infection route.

The overall prevalence of PDR among the newly confirmed HIV infections was 8.4% (303/3596; Table 2). Resistance to NNRTIs accounted for 4.7% (170/3596), which was significantly higher than resistance to protease inhibitors (PIs; 2.8%, p < 0.001) and nucleoside reverse transcriptase inhibitors (NRTIs; 1.4%, p < 0.001). Among NNRTI agents, the highest proportion of PDR was observed for nevirapine (NVP; 3.5%, 127/3596) and efavirenz (EFV; 3.4%, 121/3596), followed by RPV (2.8%, 101/3596), DOR (1.1%, 38/3596), and ETR (1.1%, 38/3596). For NRTI agents, the proportion of PDR for ABC, AZT, D4T, DDI, FTC, 3TC, and TDF was 0.7% (24/3596), 0.5% (17/3596), 0.5% (17/3596), 0.8% (30/3596), 0.4% (14/3596), 0.6% (21/3596), 0.6% (21/3596), and 0.2% (8/3596), respectively. Among PI agents, the proportion of PDR was mainly to TPV/r (1.8%, 65/3598) and NFV (0.9%, 33/3596), with PDR to other PI agents being low, around 0.1%. Notably, resistance to three or more NNRTI agents was more common (1.9%, 69/3596) than resistance to three or more NRTI agents (0.7%, p < 0.001) or PI agents (0.1%, p < 0.001). The most frequent drug resistance mutations were Q58E (1.9%, 67/3596), K103N/R/S (1.7%, 62/3596), and E138A/G/K (1.5%, 55/3596).

Factors associated with HIV PDR among the newly confirmed HIV infections are listed in Table 3. Univariable and multi-variable logistic regression analyses showed a statistically significant association between the HIV subtype and PDR. After adjusting for other factors, the individuals with the CRF07_BC subtype were more likely to cluster in transmission networks than those with the CRF08_BC subtype (aOR = 0.56, 95% CI = 0.359–0.859, p = 0.008).

The overall prevalence of ADR among the individuals with ART failure was 50.6% (83/164; Table 4). Similar to PDR, ADR to NNRTIs (46.3%, 76/164) was significantly higher than that to NRTIs (31.7%, p < 0.05) and PIs (3.7%, p < 0.001). The predominant drug resistance mutations observed were M184I/V (29.9%, 49/164), K103N (21.3%, 35/164), V179E (13.4%, 22/164), and V106I/M (9.8%, 16/164). Among the 164 individuals with ART failure, eight had viral loads exceeding 1,000 copies/mL for three consecutive years. We analyzed the dynamic changes in drug resistance mutations under drug selection pressure over these 3 years. In the first year, five NRTI resistance mutations and four NNRTI resistance mutations were detected. In the second year, one additional NRTI mutation was observed, bringing the total to six, and two additional NNRTI mutations were found, also bringing the total to six. By the third year, the number of the NRTI resistance mutations remained at six, while the NNRTI mutations increased to nine, with three new mutations added. No PI resistance mutations were detected during these 3 years. Among the NRTI resistance mutations, M184V/I (58.8%, 10/17) showed the highest mutation frequency, while among the NNRTI mutations, K103N (26.3%, 5/19), V106M (15.8%, 3/19), and G190A (15.8%, 3/19) were the most frequent. These results suggest that with prolonged antiretroviral treatment, drug resistance mutations accumulate over time, ultimately leading to clinical drug resistance.

To further investigate the characteristics of HIV drug resistance transmission, we constructed a molecular transmission network using links with a genetic distance <0.5% (Figure 1). A total of 861 individuals (23.9%, 861/3596) with 1,372 links were clustered in 278 clusters, ranging in size from 2 to 35, with 103 clusters containing three or more individuals. Among these, 54 infections with PDR were observed in 29 clusters, 10 of which included three or more individuals. We defined clusters with ≥3 nodes as large clusters.

The individuals with PDR were predominantly distributed among those infected with the CRF07_BC subtype (61.1%, 33/54), which was significantly higher than those infected with CRF01_AE (20.4%, 11/54), CRF08_BC (14.8%, 8/54), CRF55_01B (0, 0/54), URF (CRF07_BC/CRF01_AE; 1.9%, 1/54), and subtype B (1.9%, 1/54). Notably, the individuals with PDR distributed among those with CRF07_BC were more frequently found in the large clusters (60.0%, 6/10). The primary drug resistance mutations within the CRF07_BC clusters were Q58E (45.5%, 15/33), followed by K103N (15.2%, 5/33), M41L (12.1%, 4/33), and K20T (6.0%, 2/33). Fifteen individuals carrying the Q58E mutation were distributed among three clusters. The transmission fitness of the Q58E mutant strain was calculated by dividing the clustering rate of the mutant strain by the clustering rate of the wild-type strain. The transmission fitness of the Q58E mutant strain was found to be 1.02 [(15/60) / (808/3,293)].

In the CRF01_AE clusters, two individuals carried the M46I mutation, while one individual had the V108I mutation and another had the Y188F mutation; all were part of the large clusters. The individuals infected with the CRF08_BC subtype were generally older, with an average age of 51 years, and 87.5% were infected through heterosexual contact. Drug resistance mutations in the CRF08_BC clusters, including Q58E, K101E, K103N, and E138A, were all found in smaller clusters with two individuals each.

To assess the impact of ART on HIV transmission in Hangzhou, we expanded the molecular network to include both newly confirmed HIV infections and infections with ART failure (Figure 2). A total of 11 individuals with ART failure were clustered in the network. Nine of these individuals were linked to newly confirmed infections within eight clusters, and two individuals were connected to each other in one additional cluster. Among those with ART failure, 75.0% had a viral load exceeding 10,000 copies/mL compared to 25.0% with a viral load below 10,000 copies/mL. Four clusters were associated with ART failure, including one cluster that contained both newly confirmed infections and individuals with ART failure, all of whom carried the K103N drug resistance mutation.

Through the above HIV molecular transmission network analysis, 15 individuals carrying Q58E were identified within the drug-resistant transmission clusters of the CRF07_BC subtype. The clusters were named Q58E-C1, Q58E-C2, and Q58E-C3. The majority in the Q58E-C1 cluster were students (88.9%, 8/9), with an average age of 19. The majority of the individuals in the Q58E-C2 cluster lived in the Gongshu district (80.0%, 4/5). All individuals in these three clusters were men who have sex with men (MSM).

To investigate the origins of the Q58E-C1, Q58E-C2, and Q58E-C3 drug-resistant transmission clusters, Bayesian analysis was conducted to reconstruct the epidemic history of these clusters (Figure 3). We used TempEst v1.5.3 to test the molecular clock hypothesis of the selected 59 sequences. The results showed that the calculated R2 was 0.6272. The GTR model, a Log normal relaxed clock model, and a Bayesian skyline demographic growth model were used to estimate the origin time of the most recent common ancestor of the drug-resistant transmission clusters. The Skyline plot result revealed that the drug-resistant transmission clusters had an exponential growth trend from 2005 to 2010, which stabilized from 2010 to 2017 and declined in 2018. The average evolution rate of the clusters was 9.71 × 10–4 nucleotide substitutions/site/year (95% HPD: 6.70 × 10–4-1.28 × 10–3). With the support of high posterior probability, the Q58E-C1, Q58E-C2, and Q58E-C3 clusters were all introduced from Shenzhen, among which Q58E-C2 and Q58E-C3 were introduced around 2010 and showed a stable local epidemic, while the Q58E-C1 cluster was introduced around 2005 and subsequently transmitted among Hangzhou, Jiangsu, and Ningbo.