The segmented tagged-amplicon technology for HIV-1 drug resistance detection was based on the next-generation sequencing platform, and the short sequencing length was the main limitation of the NGS platform. To obtain amplicons suitable for NGS, the HIV-1 pol gene region covering the PR and first 226 codons of the RT was amplified in triple fragments with one-step RT-PCR (Takara, RR055A) and nested PCR (Takara, RR901A). The primer sequences used for amplification are shown in Supplementary Table 1. And the effective sequencing length was extended through paired-end sequencing. Each amplicon was approximately 400 bases in length, denoted pol-A, pol-B, and pol-C (nucleotide 2,169~2,583, 2,540~2,944, and 2,834~3,228 by using HXB2 as calibrator). There was a small overlap between the beginning and the end of these three amplicons. There is a small overlap between the start and end of these three amplicons, resulting in an overall amplification length of 1,060 bp. In addition, the barcode sequence was attached to the amplicons to tag different samples (Figure 1). P5/P7 tails were added at both ends of the amplified fragment to construct a high-throughput sequencing library. In a flowcell, amplicon fragment generated clusters by bridge PCR, and hundreds to thousands of amplified fragments were sequenced simultaneously. The amplicon fragments were sequenced on the Illumina sequencing platform using sequencing-by-synthesis sequencing.

In order to evaluate and validate the HIV-1 drug resistance detection method, we developed a reference group for this study. Viral load of the wild-type HIV-1 NL4-3 strain (accession number AF003887) and the resistant mutant strains were determined using quantitative real-time PCR and then mixed the wild-type HIV-1 NL4-3 strain and the resistant mutant strains according to the viral load concentration to obtain targeted DRM frequencies of 1%, 2%, 10%, 30%, and 100%. The 187 reference specimens consisted of a known ratio of virus mixtures, HIV-1 NL4-3 strain, Hepatitis B clinical samples, Hepatitis C clinical samples, and HIV-1 positive clinical samples (Supplementary Table 1). The 190 reference samples consisted of known ratios of virus mixtures, HIV-1 NL4-3 strains, Hepatitis B clinical samples, Hepatitis C clinical samples, healthy human plasma samples (HIV-1-/HBV-/HCV-), and HIV-1 positive clinical samples (Supplementary Table 2). Hepatitis B and C plasma samples were provided by clinical patients and these plasma samples were negative after the rapid HIV-1 test.

All newly HIV-1 diagnosed cases in Shenzhen from 2014 to 2015 were enrolled in this study. All patients were antiretroviral-naïve. Blood samples were collected by Shenzhen Center for Disease Prevention and Control. All participants signed written informed consents before sample collection and their background information was collected. This study was reviewed and approved by the ethics committees of the Beijing Institute of Microbiology and Epidemiology.

For the NGS, after the mixed original sequencing results were obtained on the Illumina sequencing platform, the sequencing data of different samples were split by barcode sequence. Quality control is the first step in processing high-throughput sequencing data. FastQC¹ was used for quality control of fastq files obtained from sequencing. Then, the fastq_1 and fastq_2 files obtained by paired-end sequencing were spliced to generate the sequencing data of each sample. Mothur (Schloss et al., 2009) was used to exclude poor-quality sequences from the sequencing data, and mafft (Katoh et al., 2019) was used to align the sequencing results with the reference sequence. HXB2 was used as a ruler to locate the target drug-resistant position. According to the mutations for drug resistance surveillance list of HIV drug resistance database of Stanford University², the proportion of the mutations related to drug-resistant at that location was counted.

The Sanger sequencing amplification process involved in this study was completed according to the following method. Roche High Pure Viral RNA Kit (REF:11858882001) was used to extract total viral RNA from plasma samples. The pol gene region (nucleotides positions 2,253–3,555, using HXB2 as the calibrator) spanning the protease gene and partial reverse transcriptase gene was amplified in the first round of nested reverse transcriptase-polymerase chain reaction (RT-PCR) using the One-Step RNA PCR Kit (Takara, RR055A). And the r-Taq kit (Takara, RR901A) was used in the additional PT-PCR. The primer sequences used for Sanger amplification are shown in Supplementary Table 3. Before sequencing, the PCR products were detected by 1% agarose gel electrophoresis and then subjected to direct DNA sequencing on an Applied Biosystems 3,730 Sequencer.

Sanger sequences were manually assembled in Contig Express software. The nucleotide where the secondary peak was at least 30% as high as the primary peak was counted as the degenerate bases using the International Union of Pure and Applied Chemistry (IUPAC) designations to obtain the fasta sequences. The fasta sequences were submitted to the HIV drug resistance database of Stanford University³ to generate a drug resistance report.

Due to the sensitivity of NGS to HIV low-abundance quasispecies populations, the setting of thresholds for deep sequencing data is important for reporting drug resistance variants (Kireev et al., 2018). At present, there is no accurate diagnostic value of minority drug resistance variants of HIV-1 in clinical practice. Studies have shown that the detection of NGS had a good linear relationship in the range of drug resistance mutation frequency from 1% to 100% (Lee et al., 2020). However, compared with the reporting threshold of 1%, the reporting threshold of 2% is not easy to introduce errors caused by PCR amplification and different high-throughput data analysis platforms (Chaillon et al., 2017), and has better robustness and repeatability (Becker et al., 2020; Parkin et al., 2020). Therefore, we chose a reporting threshold of 2% to describe resistance in HIV-infected individuals.

The NGS sequencing results of the review team are presented in Supplementary Table 4. We wanted to measure the overall sequence error. The overall sequence error was defined as the errors in the entire experimental process of viral nucleic acid extraction, PCR amplification, NGS sequencing, and data processing and analysis. The HIV-1 NL4-3 strain with a known sequence was used for this measurement. We used the segmented amplification method to extract, amplify, and sequence the HIV-1 NL4-3 strain and measured the average frequency of erroneous nucleotide calls. We used the well-established RNA extraction method (Roche High Pure Viral RNA Kit) to improve the RNA recovery and quality of the HIV-1 NL4-3 strain, as well as the tagged amplicon deep sequencing method to amplify and sequence. And measured the average frequency of erroneous nucleotide calls. The error rate of nucleotide sequences we measured using viral RNA ranged from 0.006% to 0.945%, with an average of 0.164%. The consensus sequences generated under the threshold greater than 1% were completely consistent with the SS sequences (Figure 2A). The detection of HIV-1 resistance mutations was the identification of amino acids, so we measured the amino acid sequence error rate of the segmented amplification method. The results showed (Figure 2B) that the average amino acid sequence error rate of pol-A gene region was 0.37% (range: 0.11%–0.92%), pol-B was 0.36% (range: 0.11%–0.83%), and pol-C was 0.61% (range: 0.16%–0.89%). This suggested that the analysis error rate of our HIVDR detection method at the amino acid level was less than 1%.

To determine the minimum viral load concentration that can be detected by our segmented amplification method, we used HIV-1 NI4-3 mutant strain as a single strain type sample, and wild-type HIV-1 NL4-3 strain and resistance mutant strain were mixed with a DRM frequency of 60% as a mixed strain type sample. Quantitative Real-time PCR was used to detect the viral load level of the samples. Each type of sample was diluted to a different viral load concentration, and each sample with a viral load concentration was divided into 5 evenly for viral nucleic acid extraction, amplification, and sequencing. We define a detection rate of less than 80% as a detection failure. It can be seen from the test results that our HIVDR detection method has a minimum viral load of 10² copies/ml for samples composed of a single strain or a sample of mixed strains (Supplementary Table 5).

Next, we were interested in whether the HIVDR detection method was consistent with the input resistance frequency and the detection frequency. Samples with different DRM frequencies were used for viral nucleic acid extraction, amplification, and sequencing. The drug resistance mutations were reported from the high-throughput sequencing data through the above-mentioned analysis process. The detection results at different frequencies were plotted as scatter plots. And we calculated the ratio of the detection frequency to the theoretical frequency. The results of the ratio were distributed around 1, which also showed that the HIVDR detection method had good accuracy (Figure 3). The average of the resistance frequency of three parallel repetitions was calculated as the actual detection frequency, and a fitting curve was made with the input theoretical frequency. From all the tested frequencies, whether it was Y181C or T219Y, the detection frequency and the theoretical frequency had highly relevant. of compliance within the range of 1%–100% (R² > 0.99), which proved that the segmented amplification method can accurately detect resistance mutations with a frequency of more than 1% (Figure 4).

The coefficient of variation was used to express the reproducibility between three replicates within a batch. We calculated the intra-batch coefficient of variation based on the results of three replicates on mixed samples with different DRM frequencies. The coefficients of variation of Y181C and T215 were between 0.08%–23.6% and 0.04%–38.22%, respectively, (Supplementary Table 6). The accuracy of this detection method was the worst when the drug resistance mutation frequency was 1%. As the frequency of drug resistance mutation increased, the precision performance was better. It was not surprising that NGS was more stable in the detection of high-frequency resistance mutations due to the sensitivity of NGS detection and the deviations produced by the PCR process.

To verify the specificity of the segmented amplification method, we selected Hepatitis B and C virus for detection. Hepatitis B and C samples were provided by confirmed clinical patients, and these samples were HIV-1 negative after rapid HIV testing. The Hepatitis B and C samples were amplified using the segmented amplification method, and the results were all negative. It showed that the amplification method had good specificity to the HIV-1 virus, and the possibility of amplifying other virus sequences was low (Supplementary Table 7).

In addition, we also used a well-characterized HIV-1 NL4-3 mutant (K101E + Y181C + T215Y) as input to verify the specificity of the segmented amplification method. The frequency of resistance mutation of three independent replicates was counted, and resistance mutations that only appear once will also be included, and the average of the reported frequency of resistance mutations that appeared multiple times will be taken (Figure 5). Judging from the HIV-1 NL4-3 (K101E + Y181C + T215Y) detection report, some ultra-low-frequency (<1%) false-positive results appeared in the results of the segmented amplification method. The detection specificity of the target mutation can reach 100%.

One hundred plasma samples of HIV-1 clinical patients were randomly trained and detected using the In-house method and segmented amplification method. The demographic characteristics of HIV-1 clinical samples are shown in Supplementary Table 8. The reads of deep sequencing data after quality control were mostly distributed around 6,000 sequences (Figure 6). A total of 10 subtypes were included in all samples. According to the detection results of surveillance drug resistance mutations (SDRMs), whether it was analyzed using our NGS-sequence analysis panel or HyDRA, the segmented amplification method can detect all the drug resistance mutations reported by Sanger sequencing. In addition, the NGS method reported an SDRM that had not been detected by SS, which indicated that the segmented amplification method had higher sensitivity than the traditional in-house method (Table 1).

Our NGS-sequence analysis panel and HyDRA were used to detect drug resistance mutations in NGS sequencing data, and the SDRMs with a frequency of more than 1% were counted (Figure 7). Due to the difference in algorithms, it was not surprising that the two methods have slight differences in the number of low-frequency SDRMs. The calculation results of McNemar’s test proved the accuracy of our NGS-sequence analysis panel for SDRM detection (p = 0.972).

All samples were amplified by PT-PCR using a segmented amplification method based on next-generation sequencing methods. Excluding factors such as insufficient plasma samples and failure of viral nucleic acid extraction, a total of 1,755 samples successfully completed the amplification and next-generation sequencing process. Among the samples for which sequencing results were obtained, 15 samples were excluded due to poor sequencing quality, and the remaining 1,740 samples obtained good sequencing results. The proportion of samples with all three segments of pol-A, pol-B, and pol-C positive was 99.15%.

A total of 1,740 participants with demographic data and NGS sequences. General characteristics of the study subjects are shown in Table 2. The age of the patients was mainly distributed between 20 and 50 years old. Most individuals were men (86.03%), single (54.25%), and Han Population (92.13%). The dominant transmission route was heterosexual contact (63.68%). The education level of junior high school and above accounted for 90.52%. CRF01_AE (37.13%) subtype was predominant, followed by CRF07_BC (38.22%) and CRF55_01B (9.54%).

According to the mutations for drug resistance surveillance in the Stanford Drug Resistance Database, a total of 360 out of 1,740 patients contained SDRMs at the 2% threshold (Figure 8). The NRTIs resistance mutations were the most prevalent (48.14%), followed by PIs (39.3%) and NNRTIs (12.56%). Thirty-five strains (2.01%) contain dual-class mutations (PI + NRTI in 24 cases, PI + NNRTI in 8 cases, and NRTI + NNRTI in 3 cases) and 2 cases contain triple-class variants. For NRTIs, the most common mutation was K65R and for PIs there were M46I + N88S. the most frequent mutations related to NNRTIs were G190E.