Written informed consents were obtained from all adults ≥ 18 years and children’s guardians ahead of sample collection. The current study was approved by the Medical Ethics Committee of Hebei provincial center for disease control and prevention [No. IRB(S)2020-031]. All experimental methods, study procedures, and study items were performed in accordance with approved regulations and guidelines.

In this study, we collected 1,008 blood samples from newly diagnosed HIV-1 individuals between 2018 and 2022. These participants’ baseline data, such as age, gender, nationality, infection routes, first CD4 cell count, marital status, education level, occupation, sample source, etc., were obtained from the HIV database of the National Center for AIDS/STD Control and Prevention. For geographic location, the participants were grouped into four regions according to the Hukou, including Jibei, Jizhong, Jidong, and Jinan. The participants should meet the following criteria: (a) age ≥15 years; (b) have sample collection time; (c) naïve treatment; and (d) newly diagnosed HIV-1 infection. A cross-sectional investigation on HIV-1 PDR was conducted, integrating detection of DR genotype, molecular networks, and statistical analysis.

Based on the method previously reported by us (Lu et al., 2024), HIV-1 RNA was extracted from 200 µL of blood plasma using the Roche MagNa pure total RNA kit (Qiagen, Valencia, CA, USA) and the pol gene fragment (HXB2:2068–5221) of HIV-1 RNA was amplified using the previous primers (Lu et al., 2024) and cycling conditions (Lu et al., 2017). Beijing Biomed Gene Technology Co., Ltd. (Beijing, China) carried out sequence determination using Sanger’s method. HIV-1 subtypes were preliminarily inferred using the online HIV Blast Microsoft (https://www.hiv.lanl.gov/content/sequence/BASIC_BLAST/basic_blast.html), followed by the phylogenetic analysis based on pol gene sequences using MEGA 7.0 and confirmed using the online REGA HIV-1 Subtyping Tool-Version 3.0 (http://dbpartners.stanford.edu:8080/RegaSubtyping/stanford-hiv/typingtool/).

HIV-1 pol sequences containing proteinase, reverse transcriptase, and integrase gene coding regions were submitted to the Stanford University HIV Drug Resistance Database (http://hivdb.stanford.edu/). DRMs to PIs, NRTIs, NNRTIs, and INSTIs were analyzed using the HIVDB algorithm version 9.5.1. Moreover, the resistance level was classified into four categories according to individual mutation score, namely, potential low-level resistance (10–14), low-level resistance (15–29), intermediate-level resistance (30–59), and high-level resistance (≥60). Molecular transmission networks were constructed based on the study pol sequences. Pairwise genetic distances were calculated using HYPHY2.2.4 with a Tamura-Nei 93 (TN93) model. We selected a genetic distance threshold of 0.015 substitutions/site to construct networks because this threshold is consistent with recent and rapid transmission. Molecular transmission networks were visualized using HIV-Trace-1.5.0.

Statistical analysis was implemented using SPSS 23.0 (SPSS Inc., Chicago, IL, USA). Means and frequencies were used to summarize demographic data. Differences in categorical variables were analyzed using the chi-square (χ²) test. The multivariable logistic regression was carried out in order to identify the association between factors of interest and the presence of PDR. Statistical significance was defined as p < 0.05 for all tests. Spearman’s method was utilized to analyze the epidemic trend of HIV-1 PDR with a p-value of 0.05 (χ² trend): r < 0 denotes a negative correlation, and r > 0 denotes a positive correlation.

In total, 1,008 participants were enrolled in this study. Of them, 925 HIV-1 pol sequences (HXB2:2068–5221) were obtained, and the PCR positive rate was 91.8%. As shown in Table 1 , men accounted for 97.4% (901/925) of the participants. Participants aged 25–49 accounted for 60.0% (555/925). The breakdown of participant occupation was as follows: farmer, 28.1% (260/925); housework/freelance, 26.5% (245/925); service, 21.2% (196/925); cadre/staff/teacher, 15.7% (145/925); and student, 8.5% (79/925). A total of 46.7% (432/925) were unmarried. Han ethnicity accounted for 939% (869/925). A total of 67.1% (621/925) had an educational level of high school or below, followed by college or above (32.9%, 304/925). A total of 55.0% (509/925) had CD4 cell counts > 200 cells/μL. Men who have sex with men (MSM) accounted for 93.0% (860/925). Among four sample sources, voluntary consultation and testing (VCT) was the most frequent, accounting for 62.2% (575/925), followed by visit test (23.2%, 215/925).

There were 77 (8.3%) of 925 study participants who had DRMs. Detailed DRMs are listed in Table 2 . There were the most locations of gene mutations in 2020 during the five years, including PI-, NRTI-, NNRTI and INSTI-DRMs. However, no PI-resistant mutation was found in 2018, 2019, 2021, and 2022, and NRTI-DRM was also not found in 2018 and 2022, respectively. HIV-1 gene mutations resistant to ARDs were found each year between 2018 and 2022 in the NNRTI and INSTI gene coding regions, respectively. Eleven of 77 participants presented multiple drug resistance: in 2019, one participant (0.1%, 1/925) had two-point DRMs in the NNRTI gene coding region; in 2020, two participants (0.2%, 2/925) had double-point DRMs in the PI gene coding region, one participant (0.1%, 1/925) had three-point DRMs in the NRTI gene coding region, and four participants (0.4%, 4/925) had double (two cases) or four (two cases)-point DRMs in the NNRTI gene coding region; in 2022, two participants (0.2%, 2/925) carried a double-class PDR mutation pattern (NNRTI + INSTI), and another participant carried three mutations (S68SG, L74M, and A98G), but both S68SG and L74M cannot cause drug resistance.

As shown in Tables 2 and 3 , in the NRTI and PI gene coding regions, T215I/TS mutations were the most common, accounting for 0.3% (3/925), followed by L210LW (0.2%, 2/925) and D67DN (0.2%, 2/925), associated with the low-level resistance to zidovudine (AZT). The NNRTI gene coding region contained the most DRM positions: E138EK/G was the most common, accounting for 1.6% (12/925), followed by K103N (1.4%, 13/925), G190GE/A/S (0.6%, 6/925), K101E (0.5%, 5/925), A98G (0.4%, 4/925), L100I (0.2%, 2/925), and Y188L (0.2%, 2/925), associated with the different level of resistance to doravirine (DOR), efavirenz (EFV), etravirine (ETR), nevirapine (NVP), and rilpivirine (RPV).

In the INSTI gene coding region, six mutations were identified, namely, four major mutations (P145PS, Q148QH, Y143S, and T66A) and two accessory mutations (S153SF and G163GRS/EK). Of them, the most frequent INSTI-DRMs were S153SF (0.6%, 6/925) and G163GRS/EK (0.6%, 6/925), followed by P145PS (0.2%, 2/925), Y143S (0.2%, 2/925), Q148QH (0.1%, 1/675), and T66A (0.1%, 1/675). G163GRS/EK, P145PS, Y143S, and T66A were associated with the resistance to elvitegravir (EVG) and raltegravir (RAL), and S153SF and Q148QH were mainly related to the resistance to dolutegravir (DTG), bictegravir (BIC), and cabotegravir (CAB).

It needs to be emphasized that the most frequent ARDs with a low- to high-level resistance were as follows: NVP (5.2%, 48/925), RPV (3.5%, 32/925), EFV (3.4%, 31/925), ETR (2.7%, 25/925), DOR (1.8%, 17/925), EVG (1.1%, 10/925), RAL (1.1%, 10/925), DTG (0.7%, 7/925), BIC (0.7%, 7/925), and CAB (0.7%, 7/925).

As shown in Figure 1 , the overall prevalence of HIV-1 PDR increased from 6.3% (8/128) in 2018 to 11.8% (25/211) in 2020, while a rapid decline to 6.4% (21/330) in 2022 was observed. The resistance prevalence of NNRTIs was the highest (5.2%, 48/925), followed by INSTIs (1.9%, 18/925), NRTIs (1.2%, 11/925), and PIs (0.2%, 2/925). In 2020, the resistance prevalence of NRTIs, PIs, and NNRTIs increased the highest value, respectively, and then declined significantly. However, the resistance prevalence of INSTIs increased from 2.3% (3/128) in 2018 to 2.7% (4/147) in 2019, while a significant decline to 0.9% (1/109) in 2021 was observed, and then increased sharply to 2.1% (7/330) in 2022. Statistical analysis ( Table 4 ) revealed that the resistance prevalence of NNRTIs, INSTIs, NRTIs, and PIs showed no increase or decrease in epidemic trend (p > 0.05); however, all kinds of PDRs presented a positive correlation (r > 0) except for NRTI-PDR (r < 0). Furthermore, Figure 1 also indicated that the prevalence of CRF01_AE and CRF07_BC showed a decreasing and increasing trend, respectively, which was associated with the spread of HIV-1 PDR strains.

In total, 15 kinds of HIV-1 subtypes were found in this study. Of them, 12 recombinant forms and three simple subtypes were identified. CRF01_AE (49.1%, 454/925) and CRF07_BC (32.6%, 302/925) were the most frequent subtypes. The prevalence of URFs (6.2%, 57/925) was significant higher than that of subtype B, becoming the third largest subtype. As indicated in Table 1 , 77 participants harboring the PDR were mainly distributed in CRF01_AE (44.2%, 34/77), CRF07_BC (41.6%, 32/77), CRF55_01B (6.5%, 5/77), B (3.9%, 3/77), CRF65_cpx (1.3%, 1/77), CRF68_01B (1.3%, 1/77), and URFs (1.3%, 1/77). All of the PI- and INSTI-DRMs were circulating in participants with CRF01_AE and CRF07_BC; however, NNRTI-DRMs were distributed in six subtypes ( Table 1 ). There were no statistical differences in resistance to PIs, NRTIs, NNRTIs, and INSTIs between different HIV-1 subtypes, ages, occupations, nationalities, educational levels, infection routes, and sample sources. However, the distribution of PI-, NRTI-, NNRTI- and INSTI-DRMs in initial CD4 counts and marital status showed obvious statistical differences (p < 0.05). The resistance rate will decline with the increase of initial CD4 counts, which suggests that early detection and early treatment can significantly reduce the resistance rate and increase therapy effect.

In total, 385 of 925 study sequences were detected within molecular transmission networks with a genetic distance threshold of 0.015. A total of 106 transmission clusters were constructed, and cluster size was 2 to 99 sequences ( Figure 2 ). Of 106 transmission clusters, 16 contained DRMs. A total of 30 resistant sequences were circulating in 16 networks with HIV-1 PDR, accounting for 62.5% of 77 total participants with DRMs. Among 30 participants with DRMs, NNRTI-DRMs were the most common mutations (n = 17), namely, E138EG/K (n = 6), K103N (n = 4), Y181C (n = 1), Y188L (n = 1), L100L (n = 1), A98G (n = 2), Y318YF (n = 1), and G190GE (n = 1), followed by 9 INSTI-DRMs [G163G/R/S (n = 5), S153SF (n = 2), T66A (n = 1), Y143S (n = 1)], 3 NRTI-DRMs [T215TI/S (n = 2) and L74LI (n = 1)], and 2 PI-DRMs [F53L/I54L (n = 1) and N88NT/L90LM (n = 1)].

Most (92.7%) of the 385 sequences included in clusters were infected through MSM, followed by heterosexual contact (7.3%). Furthermore, sequences circulating in molecular transmission networks were from all 11 cities of Hebei province. In particular, the largest cluster harbored 99 CRF_07 BC sequences, and 11 sequences had DRMs in this cluster. The largest cluster included sequences from nine cities, forming an extensive resistant transmission network. Within the largest cluster, DRM sites included G163 (n = 3), Y143 (n = 1), K103 (n = 4), E138 (n = 1), Y188 (n = 1), L100 (n = 1), T215 (n = 1), and F53/I54 (n = 1), and one participant from Qinhuangdao had 58 links.

Table 5 lists 12 potential risk factors associated with PDR to PIs, NRTIs, NNRTIs, and INSTIs. Of these 12 factors, subtype and sample source were clearly associated with HIV-1 PDR (p < 0.05). Compared with participants with CRF01_AE, those who had CRF07_BC had 1.79 times greater odds of PDR [odds ratio (OR) 1.79, 95% CI 1.021–3.168]. Compared to participants with VBD, those with VCT had 0.27 times greater odds of PDR (OR 0.27, 95% CI 1.021–3.168).