Contribution of Human Immunodeficiency Virus Type 1 Minority Variants to Reduced Drug Susceptibility in Patients on an Integrase Strand Transfer Inhibitor-Based Therapy

The role of HIV-1 minority variants on transmission, pathogenesis, and virologic failure to antiretroviral regimens has been explored; however, most studies of low-level HIV-1 drug-resistant variants have focused in single target regions. Here we used a novel HIV-1 genotypic assay based on deep sequencing, DEEPGEN (Gibson et al 2014 Antimicrob Agents Chemother 58∶2167) to simultaneously analyze the presence of minority variants carrying mutations associated with reduced susceptibility to protease (PR), reverse transcriptase (RT), and integrase strand transfer integrase inhibitors (INSTIs), as well as HIV-1 coreceptor tropism. gag-p2/NCp7/p1/p6/pol-PR/RT/INT and env/C2V3 PCR products were obtained from twelve heavily treatment-experienced patients experiencing virologic failure while participating in a 48-week dose-ranging study of elvitegravir (GS-US-183-0105). Deep sequencing results were compared with (i) virological response to treatment, (ii) genotyping based on population sequencing, (iii) phenotyping data using PhenoSense and VIRALARTS, and (iv) HIV-1 coreceptor tropism based on the phenotypic test VERITROP. Most patients failed the antiretroviral regimen with numerous pre-existing mutations in the PR and RT, and additionally newly acquired INSTI-resistance mutations as determined by population sequencing (mean 9.4, 5.3, and 1.4 PI- RTI-, and INSTI-resistance mutations, respectively). Interestingly, since DEEPGEN allows the accurate detection of amino acid substitutions at frequencies as low as 1% of the population, a series of additional drug resistance mutations were detected by deep sequencing (mean 2.5, 1.5, and 0.9, respectively). The presence of these low-abundance HIV-1 variants was associated with drug susceptibility, replicative fitness, and coreceptor tropism determined using sensitive phenotypic assays, enhancing the overall burden of resistance to all four antiretroviral drug classes. Further longitudinal studies based on deep sequencing tests will help to clarify (i) the potential impact of minority HIV-1 drug resistant variants in response to antiretroviral therapy and (ii) the importance of the detection of HIV minority variants in the clinical practice.


Introduction
Antiretroviral therapy based on the combination of several anti-HIV-1 drugs is the gold standard of care for HIV-1 infected individuals in most developed countries and has been credited with considerable reductions in morbidity and mortality [1,2]. Commonly prescribed antiretroviral regimens include two nucleoside reverse transcriptase inhibitors (NRTIs) in combination with a nonnucleoside reverse transcriptase inhibitor (NNRTI), a protease inhibitor (PI), or an integrase strand transfer inhibitor (INSTI) [3,4]. INSTI is the most recent class of antiretroviral drugs approved by the U.S. Food and Drug Administration (FDA) for the treatment of HIV-1 infection. Raltegravir (RAL, MK-0518, Isentress, Merck & Co., Inc.) was the first INSTI approved in 2007 [5]. Elvitegravir (EVG, JTK-303/GS-9137, Gilead Sciences) [6] was approved in 2012 in a fixed dose combination with a pharmacokinetic enhancer (cobicistat) and two nucleos(t)ide analog RT inhibitors (emtricitabine and tenofovir) for the treatment of antiretroviral-naïve HIV-infected individuals (QUAD, Stribild, Gilead Sciences) [7]. Dolutegravir (DTG, S/ GSK1349572, Tivicay, GlaxoSmithKline), a second-generation INSTI [8], is the latest INSTI available to treat HIV infection [9,10]. Unfortunately, despite the success of more potent and safe antiretroviral drugs with simpler regimens, current therapies have not been able to eradicate latent reservoirs of HIV-1 [11,12]. Moreover, the emergence [2,13] and potential transmission [14] of HIV-1 drug-resistant variants continues to be a cause of virologic failure in patients treated with combinations of antiretroviral drugs [2,15].
HIV-1, like other RNA viruses, has a high mutation rate (approximately 10 24 to 10 25 mutations per nucleotide and cycle of replication [16]), which coupled with selection and rapid turnover [17] results in the generation of swarms of mutants known as viral quasispecies (HIV-1 variants) [18,19,20]. HIV-1 is constantly evolving and adapting, exploring all potential combinations of mutations that could increase the capacity of the virus to replicate in any given environment (replicative fitness), including the generation and potential selection of strains carrying drug resistance mutations [21,22]. Therefore, all drug-resistant HIV-1 strains arise as the consequence of the introduction of one -or more-incorrect nucleotide during virus replication in the absence of proofreading mechanisms [23]. These drug-resistant variants will initially be present as minority members of the virus population, which could eventually be selected and outcompete other variants depending of their ability to replicate under drug pressure [21,22,24]. Current genotypic HIV-1 drug resistance assays based on population (Sanger) sequencing are able to detect these minority variants when their frequency reaches approximately 20% of the virus population [25,26,27,28,29]; however, early detection of drug resistant HIV-1 minority variants has been associated with the ability to predict clinical outcome [30,31,32,33,34,35]. Thus, ultrasensitive HIV-1 genotyping assays based on allele-specific polymerase chain reaction [36,37], oligonucleotide ligation [38,39], or deep sequencing [40,41,42,43], which are capable of detecting drug-resistant minority HIV-1 variants below the 20% level, may help to clarify the actual clinical relevance of these minority members of the viral population [30,39,44,45,46].
In this study we used a novel HIV-1 genotyping assay based on deep sequencing (DEEPGEN) [43] to detect and quantify low-level drug-resistant HIV-1 variants in twelve patients experiencing virologic failure while participating in a 48-week dose-ranging study of elvitegravir (GS-US-183-0105). Deep sequencing-based HIV-1 genotypes were then compared with (i) mutation profiles obtained by standard population sequencing, (ii) drug susceptibility using HIV-1 phenotypic assays, (iii) HIV-1 replicative fitness, and (iv) virological response to antiretroviral treatment.
Reverse transcription (RT)-PCR amplification of gag-p2/ NCp7/p1/p6/pol-PR/RT/IN-and env-C2V3-coding regions Plasma viral RNA was purified from pelleted virus particles by centrifuging one milliliter of plasma at 20,000 g660 min at 4uC, removing 860 ml of cell-free supernatant and resuspending the pellet in the remaining 140 ml, to finally extract viral RNA using QIAamp Viral RNA Mini kit (Qiagen; Valencia, CA). Viral RNA was reverse-transcribed using AccuScript High Fidelity Reverse Transcriptase (Stratagene Agilent; Santa Clara, CA) and the corresponding antisense external primers in 20 ml reaction mixture containing 1 mM dNTPs, 10 mM DTT and 10 units of RNAse inhibitor. The HIV-1 genomic region encoding the Gag proteins p2, p7, p1 and p6, and the protease, reverse transcriptase, and integrase enzymes was amplified as two overlapping fragments (1,657 nt and 2,002 nt corresponding to the p2-59half RT and 39half RT-INT, respectively) using a series of external and nested primers with defined cycling conditions [49]. External PCR reactions were carried out in a 50-ml mixture containing 0.2 mM dNTPs, 3 mM MgCl 2 and 2.5 units of Pfu Turbo DNA Polymerase (Stratagene). Nested PCR reactions were carried out in 50-ml mixture containing 0.2 mM dNTPs, 0.3 units of Pfu Turbo DNA Polymerase and 1.9 units of Taq Polymerase (Denville Scientific; Metuchen, NJ). A fragment corresponding to the C2V3 region (480 nt) of the surface glycoprotein (gp120) in the envelope gene was amplified using a series of external and nested primers with defined cycling conditions as previously described [50].
Construction of gag-p2/NCp7/p1/p6/pol-PR/RT/INrecombinant viruses Infectious recombinant viruses were constructed from each clinical samples in a HIV-1 NL4-3 backbone using a novel yeastbased cloning technology as described [49]. Briefly, PCR products spanning the 39 end of gag (p2/p7/p1/p6) and the entire pol gene (PR/RT/IN; p2-INT) were introduced via yeast homologous recombination into pRECnfl-TRPDp2-INT/URA3 or pRECnfl-TRPDINT/URA3 vectors, respectively, containing a near-full length HIV-1 genome with the yeast uracil biosynthesis (URA3) gene replacing the respective p2-INT HIV-1 coding sequences. Following yeast transformation, vector DNA was purified from the entire number of yeast colonies (typically 200 to .1,000 individual colonies) and used to transform Electrocomp TOP10 bacteria (Invitrogen). Plasmid DNA from the entire bacteria preparationto guarantee the continuity of the viral population that may have existed in vivo -was purified from 10 ml of bacteria culture (QIAprep Spin Miniprep Kit, Qiagen) and used to introduce the patient-derived HIV-1 sequences into a pNL4-3-hRluc vector expressing the human Renilla luciferase gene (hRluc) [51] as described [49]. Four micrograms of the resulting plasmid were transfected into HEK293T cells using GenDrill (BamaGen Bioscience; Gaithersburg, MD). Cell culture supernatant was harvested 48 hours post-transfection, clarified by centrifugation at 7006g, filtered through a 0.45 mm steriflip filter (Millipore; Billerica, MA), aliquoted, and stored at 280uC until further use. Tissue culture dose for 50% infectivity (TCID 50 ) was determined in triplicate for each serially diluted virus using the Reed and Muench method [52] and viral titers expressed as infectious units per milliliter (IU/ml).

Drug susceptibility determination using VIRALARTS
Drug susceptibility of twelve p2-INT-recombinant viruses was measured by determining the extent to which antiretroviral drugs inhibited viral replication in MT-4 cells as described [49]. Briefly, serial dilutions spanning empirically determined ranges of each drug were added in triplicate in 96-well plates in RPMI medium with Lglutamine (Cellgro; Mediatech) supplemented with 10% fetal bovine serum, 100 U of penicillin/mL, 100 mg of streptomycin/ mL, (Mediatech) and 10 mM HEPES (Sigma-Aldrich). MT-4 cells were infected with either the reference virus (HIV-1 NL4-3-hRluc ) [51] or the corresponding query virus (HIV-1 p2-INT-hRluc ) expressing human Renilla luciferase at a multiplicity of infection (MOI) of 0.005 IU/cell for one hour at 37uC, 5% CO 2 . HIV-infected MT-4 cells were then resuspended in RPMI medium and 30,000 cells were added to each well containing pre-plated antiretroviral drugs. Virus replication was quantified 72 hours post-infection by measuring renilla luciferase activity (relative light units, RLU) using the Renilla Luciferase Assay System (Promega, Madison, WI) in a multiwell plate reader (Victor V multilabel reader, PerkinElmer, Waltham, MA). Drug concentrations required to inhibit virus replication by 50% (EC 50 ) were calculated by (i) plotting the percent inhibition of luciferase activity versus log 10 drug concentration and (ii) fitting the inhibition curves to the data using nonlinear regression analysis (GraphPad Prism v.6.0b, GraphPad Software, La Jolla, CA). Fold change (FC) resistance values were calculated by dividing the mean EC 50 of the query virus (HIV-1 p2-INT-hRluc ) by the mean EC 50 of the internal control (HIV-1 NL4-3-hRluc ) in each assay.

HIV-1 coreceptor tropism determination using VERITROP
The ability of the virus to use the chemokine receptors CCR5 and/or CXCR4 as coreceptors to enter the host cell was quantified using a novel assay based on a modified version of the a-complementation assay for HIV-envelope glycoproteinmediated fusion [53] as previously described [54]. Briefly, patientderived PCR products spanning the gp120/gp41-coding region of HIV-1 were introduced via yeast homologous recombination into the pRECnfl-LEU-DEnv(gp120-tatex2)/URA3 vector containing a near-full length HIV-1 genome with the yeast uracil biosynthesis (URA3) gene replacing the gp120/gp41 HIV-1 coding sequence. Following yeast and bacteria transformation, 2 mg of the HIVexpression vector and 2 mg of a vector expressing the a fragment of the b-galactosidase gene (pCMVa) were co-transfected into 7610 5 HEK293T (donor) cells using Lipofectamine 2000 (Invi-trogen). The target cells (U87.CD4.CCR5 or U87.CD4.CXCR4) were transfected with 4 mg of a vector expressing the omega fragment (pCMVv) of the b-galactosidase gene. Forty-eight hours post-transfection the donor and target cells were washed, removed from the cell-culture plates, counted and re-suspended in DMEM at a concentration of 2610 6 cells per milliliter. Fifty microliters (1610 5 ) of donor and target cells were mixed and added together into 96-well plate and incubated for 4 hours at 37uC in 5% CO 2 . Cell-to-cell fusion events were quantified by measuring luminescence related to b-galactosidase activity (relative light units, RLU) using Galacto-star system (Applied Biosystems, Bedford, MA) in a multi-well plate reader (Victor V multilabel reader, PerkinElmer, Waltham, MA). Controls were run in each test, including mock cell and transfections with plasmid DNA mixtures containing (i) 100%+0%, (ii) 1%+99%, (iii) 0.3%+99.7%, and (iv) 0%+100% of vectors expressing the env gene from the 64 HIV-1 NL4-3 or the R5 HIV-1 BaL strains, respectively. Technical cutoffs for the quantification of env-mediated cell fusion events were calculated as the mean plus two standard deviations (SD) of the b-galactosidase activity detected after HEK293T cells, transfected with 100% R5 HIV-1 BaL or 100%64 HIV-1 NL4-3 , were incubated in cell-to-cell fusion experiments with U87.CD4.CXCR4 or U87.CD4.CCR5.cells, respectively.

HIV-1 replicative fitness determination using viral growth kinetics analysis
The ability of the twelve p2-INT-recombinant viruses, plus the HIV-1 NL4-3 wild-type control, to replicate in the absence of drug pressure was determined by measuring viral growth kinetics as described [49,55]. Briefly, 3610 6 MT-4 cells were infected at a MOI of 0.01 IU/cell in one ml of culture medium, incubated for 2 hrs at 37uC in 5% CO 2 . HIV-infected cells were then washed two times with 16PBS, then split to be cultured in triplicate wells of a 24-well plate (1610 6 cells/well). Culture supernatant was assayed for up to 30 days post-infection as described [56]. Viral replication was quantified using the slope of the growth curves and performing linear regression analysis derived from the equation log(y) = mt+log(h), where y is virus quantity (cpm), t is time in days, and h is the y-intercept (day 0). All slope values for each virus were used to calculate the mean, standard deviation, and 10 th & 90 th percentiles. Differences in the mean values were evaluated using a One Way Analysis of Variance test and the significance difference from the reference HIV-1 NL4-3 virus calculated using the Bonferroni's Multiple Comparison Test (GraphPad Prism v.6.0b, GraphPad Software).
Deep sequencing of gag-p2/NCp7/p1/p6/pol-PR/RT/INand env-C2V3-coding regions The three PCR products corresponding to the gag-p2/NCp7/ p1/p6/pol-PR/RT/IN-(1,657 nt and 2,002 nt fragments) and env-C2V3-(480 nt fragment) coding regions of HIV-1 were sequenced using DEEPGEN, a novel HIV-1 genotyping and coreceptor tropism assay [43]. Briefly, the three amplicons were purified (Agencourt AMPure XP, Beckman Coulter) and quantified (2100 Bioanalyzer DNA 7500, Agilent Technologies) prior to using the Ion Xpress Fragment Library Kit (Life Technologies, Carlsbad CA) to construct a multiplexed library for shotgun sequencing on the Ion Personal Genome Machine (PGM, Life Technologies). For that, a mixture of all three purified DNA amplicons (33 ng each) was randomly fragmented and blunt-ends repaired using the Ion Shear Plus Reagent (Life Technologies) followed by DNA purification (Agencourt AMPure XP, Beckman Coulter). The P1 adapter and one of 12 barcodes were ligated to the repaired fragment ends prior to DNA purification (Agencourt AMPure XP, Beckman Coulter). DNA fragments were then selected by size (i.e., 280 to 320 bp; Pippin Prep, Life Technologies) and each barcoded library, i.e., a mixture of all three amplicons per sample, was purified (Agencourt AMPure XP, Beckman Coulter) and normalized using the Ion Library Equalizer Kit (Life Technologies). All thirteen barcoded DNA libraries, corresponding to twelve patient-derived amplicons plus the HIV-

Read mapping, variant calling, and phylogenetic analysis
Reads were mapped and aligned against sample-specific reference sequences constructed for the two genomic regions (gag-p2/NCp7/p1/p6/pol-PR/RT/IN-and env-C2V3) using the DEEPGEN Software Tool Suite as described [43]. The frequency of each amino acid present in each HIV-1 genomic position was calculated and summarized in a graphical interface with particular focus on sites of known drug resistance based on the latest edition of the IAS-USA HIV Drug Resistance Mutations list [57]. A list of the amino acids at these positions, and their frequencies, was exported as a tabulated text file and used with the HIVdb Program Genotypic Resistance Interpretation Algorithm from the Stanford University HIV Drug Resistance Database (http://hivdb.stanford. edu) to infer the levels of susceptibility to protease, reverse transcriptase, and integrase inhibitors. In addition, for each dataset, reads spanning amino acid positions (i) 50 to 85 in the protease (HXB2 2,400 to 2,508), (ii) 180 to 215 in the RT (HXB2 3,087 to 3,195), (iii) 130 to 165 in the integrase (HXB2 4,617 to 4,725), and (iv) 1 to 35 in the V3 region (HXB2 7,110 to 7,217) were extracted, truncated and translated for phylogenetic analysis and HIV-1 coreceptor tropism prediction as described below. Within each dataset only one representative of any identical variant was maintained, but the overall frequency stored. All variants with a frequency .1 within the population were aligned using ClustalW [58] and phylogeny reconstructed using the neighbor-joining statistical method as implemented within MEGA 5.05 [59]. In this study, minority variants were defined as amino acid substitutions detected at .1% (based on the intrinsic error rate of the system as described [43]) and ,20% of the virus population, corresponding to those mutations that cannot be determined using population sequencing [25,26,27,28,29].
Genotypic HIV-1 coreceptor tropism determination As described for DEEPGEN [43], reads spanning amino acid positions 1 to 35 in the V3 region (HXB2 7,110 to 7,217) were extracted and truncated for HIV-1 coreceptor tropism determination using Geno2Pheno [60] with a FPR of 3.5% based on optimized cutoffs for determining HIV-1 coreceptor usage as previously described [33,34,61]. Deep sequencing V3 sequences usually spanned 105 nucleotides (35 amino acids), with some minor discrepancies associated with natural HIV-1 variation [62,63], which led to V3 sequences with an open reading frame of 96, 99, 102, 108, or 111 nucleotides, all starting and ending with a cysteine codon, i.e., TG(T/C). V3 reads with stop codons (TGA, TAA, or TAG) and/or where the nucleotide length was not a multiple of 3 (e. g., 101, 103, 104, 106, etc.), mostly associated with natural or methodology (PCR or sequencing)-induced insertions and/or deletions, were not included in the analysis. Plasma samples were classified as containing non-R5 viruses if at least 2% of the individual sequences, as determined by deep sequencing, were predicted to be non-R5 [33,34].

Statistical Analyses
Descriptive results are expressed as median values, standard deviations, and confidence intervals. Pearson's correlation coefficient was used to determine the strength of association between categorical variables. All differences with a P value of ,0.05 were considered statistically significant. As described above, differences in the mean of the slope values for the viral growth kinetics curves were determined using a One Way Analysis of Variance test and the significance difference from the reference HIV-1 NL4-3 virus calculated using the Bonferroni's Multiple Comparison Test. All statistical analyses were performed using GraphPad Prism v.6.0b (GraphPad Software, La Jolla, CA) unless otherwise specified. gag-p2/NCp7/p1/p6/pol-PR/RT/IN and/or env-V3 nucleotide sequences obtained by deep sequencing in this study have been submitted to the Los Alamos National Laboratory HIV-DB Next Generation Sequence Archive (http://www.hiv.lanl.gov/content/ sequence/HIV/NextGenArchive/Gibson2014).

Results
Antiretroviral drug susceptibility of multidrug-resistant viruses determined using HIV-1 genotyping (based on Sanger sequencing) and phenotyping assays Plasma samples from 12 HIV-infected individuals experiencing virologic failure while participating in a 48-week dose-ranging study of elvitegravir were originally analyzed using standard (Sanger) population sequencing and an HIV-1 phenotyping assay (PhenoSense, Monogram Biosciences) [47,48]. Here, we used patient-derived PCR products from the same plasma samples to construct gag-p2/NCp7/p1/p6/pol-PR/RT/IN-recombinant viruses and assess their susceptibility to 22 antiretroviral drugs, including EVG and RAL, using an alternative HIV-1 phenotyping assay (VIRALARTS) [49] (Table 1). A complete list of all amino acid substitutions in the protease, RT, and integrase coding regions, determined directly from plasma-derived amplicons using Sanger sequencing, is included in Table S1. Most of these highly treatment-experienced patients failed the antiretroviral regimen with viruses carrying multiple drug resistance mutations to protease (mean 4; range 0 to 10), nucleoside reverse transcriptase (3.8, 0 to 7), nonnucleoside reverse transcriptase (0.8, 0 to 3), and integrase (2.3, 1 to 4) inhibitors ( Table 1). Most of the resistance mutations in PR and RT were present prior to initiating therapy with the EVG-based regimen (data not shown). The reduction in susceptibility to several antiretroviral drugs was corroborated by both HIV-1 phenotypic assays. For that, we compared the mutation scoring generated by the HIVdb Program Genotypic Resistance Interpretation Algorithm from the Stanford University HIV Drug Resistance Database (http://hivdb.stanford.edu) for each drug with the fold-changes in EC 50 values (FC) obtained with the HIV-1 phenotyping assays. In general, there was good agreement between the Sanger sequencing-based genotype and the PhenoSense or VIRALARTS tests (r values of 0.89 or 0.87, respectively, p,0.0001 Pearson coefficient correlation). As expected, a strong statistically significant correlation was observed between the EC 50 fold changes obtained with both HIV-1 phenotyping assays (r = 0.99, p,0.0001 Pearson coefficient correlation).
Replicative fitness of multidrug-resistant gag-p2/NCp7/ p1/p6/pol-PR/RT/IN-recombinant viruses All 12 replication-competent gag-p2/NCp7/p1/p6/pol-PR/ RT/IN-recombinant viruses, generated for drug susceptibility analysis with VIRALARTS, were used in viral growth kinetic experiments to determine their replicative fitness in the absence of drug pressure. Statistical analysis of the slope of the growth curves showed that the replicative fitness of most multidrug-resistant viruses was statistically significantly reduced compared to the reference HIV-1 NL4-3 strain (Fig. 1). The differences in replicative fitness among the recombinant viruses were patient-dependent and guided mainly by the number and type of drug resistance mutations in the PR, RT, and INT coding regions ( Fig. 1B and Table S1). Interestingly, the replicative fitness of the 08-198 recombinant virus was not impaired relative to the wild type HIV-1 NL4-3 strain, despite carrying primary drug-resistance mutations in the reverse transcriptase (M184V + T215F) and integrase (T66A + E92Q + S147G) coding regions (Fig. 1B). Nevertheless, there was a tendency -although not statistically significant-for viruses with higher mutation scoring determined by Sanger sequencing to have reduced viral replicative fitness values in the absence of drug pressure (r = 20.2, p = 0.16, Pearson coefficient correlation). Similar results were observed when viral replicative fitness values were compared with EC 50 fold changes determined with either one of the HIV-1 phenotyping assays, i.e., PhenoSense or VIRALARTS (r values of 20.14, p = 0.23 or 20.13, p = 0.21, respectively, Pearson coefficient correlation).

Determination of antiretroviral drug susceptibility using deep sequencing
We have recently developed a novel HIV-1 genotyping and coreceptor assay based on deep sequencing (DEEPGEN) that allows the detection of minority HIV-1 variants when present at frequencies as low as 1% of the HIV-1 population [43]. Here we used this methodology to identify low-level drug-resistant viruses otherwise not detected by Sanger sequencing. All 12 samples were multiplexed into a single Ion 316 chip (59% loading efficiency of Ion Sphere Particles), generating 2,581,962 total quality reads and an average read length of 176 bp. Although comparable, the average sequencing coverage at each nucleotide position varied with each sample and HIV-1 genomic region analyzed, i.e., gag-p2/NCp7/p1/p6 (mean 8,940; range 3,881 to 11,647), protease (10,389; 5,319 to 13,299), reverse transcriptase (9,426; 4,389 to 18,055), integrase (6,562; 2,512 to 8,894), and env-V3 region (2,310; 1,026 to 3,270) (Fig. 2). These metrics ensured the minimum coverage of 1,000 per nucleotide position sequenced required to guarantee the detection of a minor variant present at least at 1% of the population [64].
We first constructed neighbor-joining phylogenetic trees using 105 bp fragments from the protease, RT, integrase, or V3 HIV-1 regions -obtained by Sanger sequencing-to rule out any potential cross-contamination (Fig. 3A). Next, all deep sequencing reads with a frequency .1 corresponding to the same protease, RT, integrase, and V3 HIV-1 regions were used to construct neighborjoining phylogenetic trees to quantify intra-and inter-patient genetic distances. A total of 2,482 unique protease (mean 207; range 80 to 581), 2,325 unique reverse transcriptase (194; 64 to 730), 1,462 unique integrase (122; 64 to 186), and 695 unique V3 (58; 22 to 75) sequences were included in each phylogenetic analysis. A clear virus-dependent clustering of unique sequences was evident only in the V3 region, while certain protease, RT, and integrase HIV-1 sequences from different patients branched together due to the multiple drug resistance mutations shared by the viruses obtained from these highly antiretroviral-experienced patients (Fig. 3B). Nevertheless, and as expected, interpatient genetic distances were larger than the range of intrapatient genetic diversity in the four HIV-1 regions, i.e., 0.0662 (0.0082 to 0.0264), 0.1255 (0.0106 to 0.0418), 0.0313 (0.0105 to 0.0276), and 0.1752 (0.0235 to 0.0738) substitutions per nucleotide in the protease, RT, integrase, and V3, respectively (Fig. 3C). Although not statistically significant, there was a tendency for viruses with higher genetic diversity in all four HIV-1 coding regions analyzed, i.e., more heterogeneous virus population, to have lower viral replicative fitness values (regression values ranging from -0.16 to -0.24, p. 0.1, Pearson coefficient correlation). Interestingly, plasma viral load seems to be related to HIV-1 genetic diversity in the protease, RT, and integrase (r = 0.26, 0.26, and 0.18, respectively, p.0.1, Pearson coefficient correlation) but may correlate negatively with quasispecies heterogeneity as measured in the V3 region (r = 20.23, p = 0.19, Pearson coefficient correlation).
Altogether a total of 194 mutations in positions associated with drug resistance were detected by both Sanger and deep sequencing, i.e., 113 in the protease, 64 in the RT, and 17 in the integrase (Fig. 4). As expected, all the drug resistance mutations identified by population sequencing were also detected by deep sequencing, while 59 additional drug resistance mutations were detected only by deep sequencing, i.e., 30 in the protease, 18 in the RT, and 11 in the integrase ( Fig. 4 and Fig. S1). Overall, the difference in the numbers of drug resistance mutations detected by both methods was significant, even when the mutations were quantified by drug class, i.e., an average of 2.5, 1.5, and 0.9 additional mutations associated with PI, RTI, and INI, respectively, were detected by deep sequencing compared to population sequencing (Paired t test, p,0.0001) (Fig. 4). Figure 5 shows a comparison of the frequency of amino acids detected by population and deep sequencing in a virus carrying multiple mutations associated with resistance to PI, RTI, and INI (08-189). A number of amino acid substitutions associated with reduced susceptibility to antiretroviral drugs were identified only by DEEPGEN, all of them at frequencies below the limit of detection (,20%) of population sequencing [25,26,27,28,29], e.g., V82A (3.5%) in the protease, D67N (16.4%) in the RT, and E92Q (4%) in the integrase coding region. As exemplified in the integrase coding region, deep sequencing not only detected additional INSTI-resistance mutations in virus populations from 6 of the 12 HIV-infected individuals but was also able to quantify mixtures of amino acid substitutions identified qualitatively by population sequencing, e.g., E92Q (30.3%), N155H (59.9%) versus E92E/Q and N155N/H in patient 08-202, respectively ( Table 2).
As described above, we used the mutation scoring generated by the HIVdb Program Genotypic Resistance Interpretation Algorithm to compare the results obtained by deep sequencing with standard HIV-1 genotyping and phenotyping data, as well as clinical parameters. As expected, a strong statistically significant correlation was observed between the Sanger-and deep sequencing-based genotypes (r = 0.95, p,0.0001 Pearson coefficient

HIV-1 coreceptor tropism determination
In addition to identification of low-level drug resistance mutations, DEEPGEN is able to detect minority non-R5 (CXCR4-, dual-tropic, and/or dual mixed) HIV-1 variants [43]. Thus, here we used this deep sequencing-based assay to determine the coreceptor tropism of viruses obtained from the 12 HIVinfected individuals and compared the results with a phenotypic HIV-1 tropism test (VERITROP) [54]. Hierarchical clustering analysis grouped the two HIV-1 coreceptor tropism determinations based on their ability to detect R5 and non-R5 sequences, with 6 out of 12 patients harboring R5 HIV-1 strains based on deep sequencing analysis of the V3 region (Fig. 6A). In these patients, the concordance was high (81.8%) between the two HIV-1 tropism methods (Fig. 6B). Although not statistically significant (Paired t test p values ranging from 0.27 to 0.41), there was a tendency for non-R5 viruses to have a higher overall drug resistance level, based on the sum of the mutation scoring for all antiretroviral drugs (Sanger or deep sequencing) or EC 50 foldchange values (PhenoSense OR VIRALARTS). Interestingly, the opposite was observed when HIV-1 coreceptor tropism was compared to plasma viral load, i.e., patients infected with R5 viruses seemed to have higher (while not significantly different), viral loads than patients carrying non-R5 viruses as majority members of the population (mean 359,150 copies/ml vs. 177,433 copies/ml, respectively; Paired t test, p = 0.101).

Discussion
The clinical significance of minority HIV-1 drug-resistant variants is still under debate [30,32,39,44,45,46,65]. While a few studies have failed to establish an association between the presence of minority variants carrying drug resistance mutations and therapeutic failure [66,67,68,69], many others have associated the identification of low-abundance drug-resistant viral variants as having an impact in treatment outcome [32,35,44,45,65,70]. It is possible that the threshold for the identification of significant low-level variants may be mutation and/or antiretroviral drug dependent [31]; thus, it is clear that additional studies based on reliable ultrasensitive assays are needed to better understand the clinical relevance of these minority HIV-1 variants. Here we used a novel HIV-1 genotyping and coreceptor tropism assay based on deep sequencing [43] to quantify minority HIV-1 variants in patients experiencing virologic failure while participating in a 48week dose-ranging study of elvitegravir and evaluated their Neighbor-joining phylogenetic trees constructed using population (Sanger) sequencing of 105-bp fragments corresponding to the HIV-1 protease, RT, integrase, and V3 regions from the 12 patients. Phylogenetic trees were rooted using the HIV-1 HXB2 sequence (GenBank accession number AF033819). (B) Neighbor-joining phylogenetic trees constructed using reads with a frequency .1 corresponding to 105-bp fragments from the protease, RT, integrase, and V3 regions. Each color-coded dot represents a unique variant, frequency is not depicted. Bootstrap resampling (1,000 data sets) of the multiple alignments tested the statistical robustness of the trees, with percentage values above 75% indicated by an asterisk. s/nt, substitutions per nucleotide. (C) HIV-1 intra-and inter-patient genetic diversities were determined using MEGA 5.05 [59]. doi:10.1371/journal.pone.0104512.g003 contribution to HIV-1 replicative fitness and virological response to the antiretroviral treatment.
All 12 treatment-experienced HIV-infected individuals failed the INSTI-based therapy with viruses carrying a multitude of mutations conferring resistance to many PIs, NRTIs, NNRTI, and/or INSTIs. Among the primary INSTI-resistance mutations, standard population sequencing identified E92Q, N155H, and Q148R in 6, 5, and 2 patients, respectively. This is in agreement with previous studies showing that E92Q, Q148R/H/K, and N155H have been the most common EVG-resistance mutations identified in vivo [55,71]. As expected, and similar to prior studies [49,72], there was a good correlation between HIV-1 drug resistance determined by genotyping (Sanger sequencing) and phenotyping (PhenoSense or VIRALARTS) assays. Nonetheless, using a deep sequencing-based HIV-1 genotyping assay (DEEP-GEN) [43] we were able to identify low-level HIV-1 variants (at a frequency of ,20% of the population) carrying 59 additional drug-resistance mutations in the protease, reverse transcriptase and integrase coding regions of viruses from these highly antiretroviral-experienced patients. Previous studies have de-  scribed similar findings, although analyzing no more than one or two drug-targeted regions at a time, in antiretroviral-naïve [65,66,73,74,75,76] or -experienced HIV-infected individuals [43,68,69,77,78,79,80], with the detection of minority mutations usually correlating with historical antiretroviral treatment. In the case of INSTI-resistance mutations, using deep sequencing we corroborated and quantified all the amino acid substitutions identified by Sanger sequencing, even in positions labeled as mixtures, e.g., E92E/Q (55.1% E92Q) or N155N/H (79.9% N155H) in patients 08-172 and 08-177, respectively. More importantly, we identified 9 additional INSTI-resistance mutations at frequencies between 2.7% and 11.6% in 6 of the 12 patients, which most likely were or would have been playing a role in the overall susceptibility to INSTIs if the patients had continued in the antiretroviral regimen. Are low-abundance HIV-1 drug-resistant variants contributing to virologic failure? As described above, additional minority drugresistance mutations detected by deep sequencing have been shown to correlate with treatment failure [46,70,75,76,77]. Here we evaluated the relationship between the presence of drugresistant variants, identified by Sanger or deep sequencing, and different virological and clinical parameters. First, the mutation scoring generated by the HIVdb Program using Sanger sequencing correlated slightly better with PhenoSense than with the other HIV-1 phenotyping test (VIRALARTS). On the other hand, a better agreement was observed between deep sequencing and VIRALARTS. VIRALARTS uses a yeast-based cloning methodology, which provides a better representation of the HIV-1 intrapatient population that other cloning methodologies [49,81]. Thus, it is possible that additional minority HIV-1 variants, detected only by deep sequencing, were introduced during the yeast-based cloning step to generate the HIV-1 recombinant viruses. Second, although not statistically significant, viruses with higher drug-resistance score (determined by deep sequencing) had lower replicative fitness values, as it is usually the case for HIV-1 drug-resistant strains [21,22,24]. Interestingly, plasma viral load correlated negatively with the level of HIV-1 drug resistance determined not only by deep sequencing but also by the other three HIV-1 genotyping or phenotyping methods (r = 20.62 to 273, p,0.0001 Pearson coefficient correlation). The accumulation of primary drug-resistance mutations, although increases the replicative fitness of the virus in the presence of drug pressure, it does affect viral replication in the absence of antiretroviral drugs [21,49,82,83]. Thus, although all patients experienced virologic failure with high plasma viral loads (mean, 268,291; range 35,400 tot 578,000 copies/ml), the addition of drug-resistance mutations seems to have a negative effect in overall viral replication. This finding was supported by the fact that plasma viral load correlated positively with HIV-1 genetic diversity in the protease, RT, and integrase (driven by the accumulation of drug resistance mutations) but it seems to correlate negatively with the env-V3 quasispecies heterogeneity, perhaps due to a bottleneck effect during the selection of HIV-1 drug-resistant variants [84,85].
In summary, here we used a novel HIV-1 genotyping and coreceptor assay based on deep sequencing to analyze the contribution of minority HIV-1 drug-resistant variants in patients experiencing virologic failure while participating in a 48-week dose-ranging study of elvitegravir. We showed that these lowfrequency drug-resistant viruses could enhance the overall burden of resistance not only to INSTI but also to PI, NRTI, and NNRTI. Should minority HIV-1 drug-resistant variants be considered when planning first-line or subsequent antiretroviral regimens? Detection of low-abundance NNRTI-or NRTI-resistant HIV-1 variants prior the initiation of antiretroviral therapy seems to correlate with a higher risk of virologic failure [32,44,45,65,70]. On the other hand, similar studies have not been able to associate the presence at baseline of low-frequency HIV-1 variants resistant to PI [67], NRTI [68,69], or NNRTI [66] with antiretroviral Figure 6. HIV-1 coreceptor tropism determination using deep sequencing (DEEPGEN [43]) and a phenotypic assay (VERITROP [54]). (A) Hierarchical clustering analysis was used to group the two HIV-1 coreceptor tropism determinations by similarity. Dendograms were calculated using the Euclidean distance and Complete cluster methods with 100 bootstrap iterations as described (http://www.hiv.lanl.gov/content/sequence/ HEATMAP/heatmap.html). Bootstrap values are indicated. Green and red blocks indicate the absence or presence of non-R5 (X4) viruses, respectively, as determined by each assay. (B) Concordance between DEEPGEN and VERITROP. doi:10.1371/journal.pone.0104512.g006 therapy failure. To the best of our knowledge, no study has used deep sequencing to associate the effect of low-level INSTIresistance variants prior the initiation of an INSTI-based antiretroviral regimen. Therefore, further and well-controlled longitudinal studies based on ultrasensitive HIV-1 genotyping assays are necessary to understand the clinical implications of minority HIV-1 drug resistance mutations.