Spectrum of Atazanavir-Selected Protease Inhibitor-Resistance Mutations

Ritonavir-boosted atazanavir is an option for second-line therapy in low- and middle-income countries (LMICs). We analyzed publicly available HIV-1 protease sequences from previously PI-naïve patients with virological failure (VF) following treatment with atazanavir. Overall, 1497 patient sequences were identified, including 740 reported in 27 published studies and 757 from datasets assembled for this analysis. A total of 63% of patients received boosted atazanavir. A total of 38% had non-subtype B viruses. A total of 264 (18%) sequences had a PI drug-resistance mutation (DRM) defined as having a Stanford HIV Drug Resistance Database mutation penalty score. Among sequences with a DRM, nine major DRMs had a prevalence >5%: I50L (34%), M46I (33%), V82A (22%), L90M (19%), I54V (16%), N88S (10%), M46L (8%), V32I (6%), and I84V (6%). Common accessory DRMs were L33F (21%), Q58E (16%), K20T (14%), G73S (12%), L10F (10%), F53L (10%), K43T (9%), and L24I (6%). A novel nonpolymorphic mutation, L89T occurred in 8.4% of non-subtype B, but in only 0.4% of subtype B sequences. The 264 sequences included 3 (1.1%) interpreted as causing high-level, 14 (5.3%) as causing intermediate, and 27 (10.2%) as causing low-level darunavir resistance. Atazanavir selects for nine major and eight accessory DRMs, and one novel nonpolymorphic mutation occurring primarily in non-B sequences. Atazanavir-selected mutations confer low-levels of darunavir cross resistance. Clinical studies, however, are required to determine the optimal boosted PI to use for second-line and potentially later line therapy in LMICs.


Introduction
Ritonavir-boosted atazanavir has become increasingly important as an option for second-line therapy in low-and middle-income countries (LMICs) [1]. Although it appears to have comparable efficacy to ritonavir-boosted lopinavir (lopinavir/r) [2,3], there are few data on the mutations arising in patients receiving boosted or unboosted atazanavir . Approximately 99% of sequences were obtained from plasma and 1% from peripheral blood mononuclear cells (PBMCs). Next-generation sequencing (NGS) was performed in 1 of the 30 studies. The most common subtypes were B (61.9%), C (13.6%), A (6.7%), G (6.1%), 02_AG (4.9%), F (3.1%) and D (1.1%). Of 1497 patients, 62.7% (n = 939) received boosted atazanavir and 37.3% (n = 558) received unboosted atazanavir. A higher proportion of patients with subtype B (70.4% of 927) compared with non-subtype B (50.2% of 570) viruses received boosted (p < 0.001). Table 2 summarizes the numbers of patients according to the administration of atazanavir (boosted vs. unboosted), subtype (B vs. non-subtype B), previous antiretroviral therapy (ART) (naïve vs. experienced), and year of ART initiation.  Footnotes: 1 DRMs were defined as those with an HIVDB drug resistance program penalty score for ≥1 PI. 2 Patients with available year of ART initiation (n = 1127) were grouped into four time periods containing approximately equal numbers of patients.
Sequences from patients receiving unboosted atazanavir were also slightly more likely to have one or more NP-TSMs compared with sequences from patients receiving boosted atazanavir (12.2% of 558 vs. 8.6% of 939; p = 0.03). Each of the 33 reported NP-TSMs occurred in similar proportions in patients receiving unboosted and boosted atazanavir. An additional 197 mutations, previously classified as nonpolymorphic treatment selected mutations (NP-TSMs), occurred in 149 sequences, including in 109 of the 264 sequences containing a PI-associated DRM and 40 of the 1215 sequences without a PI-associated DRM. There were 33 different NP-TSMs of which the most common were L89T (34.9% of 149 sequences), K55R (15.4%), I85V (11.4%), A71I (9.4%), and E34Q (7.4%) ( Table   50L 46I

Subtypes
The proportion of sequences containing one or more PI-associated DRMs was similar in subtype B (20.0% of 570) versus non-subtype B (16.2% of 927; p = 0.07) sequences (Table 2). Of the 48 reported PI-associated DRMs, G73S was significantly more common in subtype B (3.1% of 927) than non-subtype B (0.4% of 570; adjusted p = 0.005) sequences. Of the 33 reported NP-TSMs, only L/M89T was significantly more common in non-subtype B than in subtype B sequences (8.4% of 570 vs. 0.4% of 927; adjusted p < 0.001). In subtypes A, C, G, CRF01_AE, and CRF02_AG, the consensus amino acid at position 89 is methionine (M) [41] and 89T requires just a single transition in these subtypes (ATG => ACG). In contrast, changing to 89T requires a one transition plus one transversion change in subtype B (CTN or TTR => ACN).
Among the 590 previously ART-experienced patients, atazanavir was administered with 2 NRTIs in 345 (58.5%) patients. Among the remaining 245 patients, the co-administered ARVs were not provided for 163 (27.6%), while 82 (13.9%) received a variety of other ARVs. Only four patients received atazanavir plus one additional ARV.
The year of ART-initiation was available for 1127 (75.3%) of all patients. The patients could be pooled into four time periods containing approximately equal numbers spanning the years between 1993 and 2018 ( Table 2). The proportion of patients with one or more PIassociated DRMs decreased over time (binomial coefficient = −0.26; 95% CI: −0.45 to −0.07; p = 0.007), but the number of DRMs in patients with one or more DRMs did not change.
Using just those patients for whom the year of ART initiation was available, a multivariate logistic regression analysis was performed to assess the association between four factors and the development of a PI-associated DRM. The four factors included the year of ART initiation, subtype (B vs. non-subtype B), the use boosted vs. unboosted atazanavir, and previous ART (naïve vs. experienced). The analysis found that a later year of ART initiation (OR: 0.62; 95%CI: 0.49-0.79; p = 0.0001) and the administration of boosted atazanavir (OR: 0.57; 95%CI: 0.35-0.93; p = 0.02) were associated with a decreased risk of developing a PI-associated DRM.

Bayesian Network Analysis of Correlated Mutations
We used the 1437 (96%) sequences containing 0 to 4 PI-associated DRMs (i.e., sequences with ≥5 PI-associated DRMs were excluded) to calculate Jaccard similarity coefficients and their standard Z scores for all pairs of DRMs and NP-TSMs. Eleven pairs of mutations comprising six major DRMs (M46I, I50L, I54V, V82A, N88S and L90M), three accessory DRMs (K20T, L33F and G73S), and the NP-TSM L89T participated in one or more significant pairwise correlations (p < 0.01). We then performed a Bayesian network analysis to determine the conditional dependency between the mutations in each of the pairwise correlations ( Figure 2). quences with ≥ 5 PI-associated DRMs were excluded) to calculate Jaccard similarity coefficients and their standard Z scores for all pairs of DRMs and NP-TSMs. Eleven pairs of mutations comprising six major DRMs (M46I, I50L, I54V, V82A, N88S and L90M), three accessory DRMs (K20T, L33F and G73S), and the NP-TSM L89T participated in one or more significant pairwise correlations (p < 0.01). We then performed a Bayesian network analysis to determine the conditional dependency between the mutations in each of the pairwise correlations ( Figure 2).  The Bayesian network analysis yielded 11 mutation pairs, including 6 major DRMs (red), 3 accessory DRMs (yellow), and an additional nonpolymorphic treatment-selected mutation (light blue) with a significant Jaccard correlation coefficient (p < 0.01). The thickness of the arrows indicates the strength of the probabilistic relationship of the two mutations. The direction of the probabilistic causation is shown with an arrowhead. For the direction between V82A and I54V for which the probabilistic causation is not greater than the probabilistic causation of the opposite direction by 0.1, the arrowhead is not shown.

Virological Failure with Resistance
Five of the thirty studies included participants from three clinical trials and from two clinical cohorts for which genotypic resistance testing was routinely available (Table 1). Together, these five studies included 1037 (69.3%) of all 1497 patients from whom sequences were available. Of these 1037 patients, 63.0% and 37.0% received boosted and unboosted atazanavir, respectively. In these studies, the proportion of sequences containing one or more PI-associated DRMs ranged from 2.9% to 14.5% and the overall proportion of sequences containing one or more PI-associated DRMs in patients receiving boosted and unboosted atazanavir were 7.2% and 13.5%, respectively.

Studies Not Included in the Analysis
We identified 32 additional studies reporting sequences from 1089 previously PInaïve patients receiving boosted or unboosted atazanavir-containing regimens (Table S3). Approximately 10% of the sequences in these studies were reported to contain one or more PI-associated DRMs. However, as the sequences were not available and as different mutations were reported in different studies, we did not include the data from these studies in our analysis.

Discussion
The spectrum of atazanavir-selected mutations has been largely influenced by data published in the earliest in vitro passage experiments and clinical trials. During in vitro passage experiments with three subtype B clones, the most commonly emerging DRMs were V32I, M46I, I50L, I84V, and N88S [42]. The initial reports of the in vivo selection of PI-associated DRMs, based on the use of unboosted atazanavir in ART-naïve patients, demonstrated that I50L and G73S were the most commonly occurring mutations in patients with VF [24,43]. A few cases of VF and emergent PI-associated DRMs have been reported in the clinical trials of ART-naïve patients receiving boosted atazanavir [16,44], consistent with the hypothesis that PI-resistance mutations develop only in viruses exposed to a narrow window of suboptimal drug concentrations that both exert selective pressure on the virus and allow virus replication [45]. Nonetheless, PI resistance in previously PInaïve patients receiving lopinavir/r for second-line therapy has increasingly been reported, usually beginning after 12-18 months of therapy [46]. In addition, phenotypic studies have shown that DRMs selected by other PIs confer atazanavir cross resistance particularly when they occur in combination [47,48].
In the years since atazanavir has been introduced, there has been a gradual accumulation of data on the spectrum of mutations emerging in previously PI-naïve patients with VF on an ART-regimen containing boosted and less commonly unboosted atazanavir. In contrast to the earliest clinical trials of boosted atazanavir, these studies included cohorts of patients who were ART experienced at the time atazanavir was administered and who may not have been monitored as closely for VF thus enabling their viruses to evolve for longer period of time under atazanavir selection pressure. Moreover, these studies have included an increasing proportion of sequences from patients with non-subtype B viruses.
Our analysis confirmed that the five major DRMs selected in vitro by atazanavir (V32I, M46I, I50L, I84V, and N88S) were among the most commonly occurring major DRMs. Four additional DRMs also occurred commonly, including M46L, I54V, V82A, and L90M. I50L is a signature atazanavir-associated DRM because it has only been reported in patients receiving atazanavir and it increases susceptibility to PIs other than atazanavir [48,49]. N88S is also considered a signature atazanavir-associated DRM because it is rarely selected by other PIs and it does not significantly reduce susceptibility to other PIs, with the exception of nelfinavir and indinavir [48]. L/M89T may also be a signature atazanavir-associated mutation because it appears to occur more commonly in patients receiving atazanavir than in patients receiving any other PI; for example, it has only been reported in three previously PI-naïve patients receiving lopinavir/r (https://hivdb.stanford.edu/cgi-bin/Mutations. cgi?Gene=PR; accessed on 1 November 2021). In contrast, each of the remaining atazanavirselected mutations appear to be commonly selected by other PIs, in particular lopinavir/r.
With the exception of L10F and L33F, each of the above 17 most commonly selected major or accessory DRMs was significantly associated with reduced atazanavir susceptibility in a previously published weighted least squares regression analysis of 1644 sequences [48]. Few published phenotypic data are available on sequences containing L89T.
Some limitations of our review should be discussed. First, most of the sequences that we reviewed were obtained from retrospective cohort studies and case series. For these studies, the duration of therapy, accompanying ARVs, frequency of virological monitoring and genotypic resistance testing, and duration of virological failure were generally not available. Therefore, the extent to which these factors were associated with emergent PI-associated DRMs could not be explored. Second, the dataset contained an underrepresentation of subtypes other than subtype B. Third, we could not be sure that every sequence was obtained from a patient receiving atazanavir as his/her first PI as treatment histories are often incomplete. Nonetheless, we emailed the authors of those studies containing the largest numbers of DRMs and received confirmation that, to the authors' knowledge, the patients had just received atazanavir. Fourth, at least 32 studies in PubMed that contained sequences from 1000 patients receiving boosted or unboosted atazanavir could not be included in our analysis because the primary sequence data and complete list of protease mutations were not available.
In conclusion, to our knowledge, this is the only comprehensive analysis of atazanavirselected mutations. Our analysis shows that the spectrum of atazanavir-selected mutations extends beyond those mutations observed in the earliest clinical trials in which patients received either boosted or unboosted atazanavir. The expanded spectrum is likely due to the large number of sequences in our analysis and the likelihood that many of the patients in the studies we reviewed had prolonged VF and ongoing replication while receiving atazanavir. The study also identified one novel nonpolymorphic atazanavirselected mutation that predominantly occurred in non-subtype B sequences. The relatively low cross-resistance to darunavir/r combined with preliminary data suggests that boosted atazanavir can be an efficacious regimen for second-line therapy. However, comparative clinical trials are required to determine the optimal boosted PI to use for second-line and potentially later-line therapy in LMICs.

Study Selection Criteria
We analyzed publicly available HIV-1 group M protease nucleotide sequences obtained from previously PI-naïve patients receiving boosted or unboosted atazanavir. Sequences were obtained from HIVDB, which is populated with sequences from GenBank annotated with the ART history of the patients from whom the sequences were obtained [13]. The analysis was last updated 31 December 2021. We supplemented the data in HIVDB with previously unpublished sequences performed at SUH and with previously unpublished sequences from two collaborating research groups: the EIDB [14] and the RHIVDB [15]. Additionally, we performed a PubMed search to identify studies describing HIV-1 group M protease sequences that were not present either in HIVDB or GenBank.
Publications reporting eligible protease sequences were reviewed to determine the treatment history of the patient from whom each sequence was obtained to confirm that the patient had received no PI prior to atazanavir and to distinguish those patients receiving unboosted atazanavir from boosted atazanavir. Each sequence was annotated with the year and country of virus isolation, the type of sample (e.g., PBMCs), the sequencing method (Sanger dideoxynucleoside sequencing versus NGS), and the nature of the study population. HIV-1 subtype was determined using the HIVDB subtyping program [50].
We also characterized each study according to whether it included patients in a clinical trial or in a treatment cohort for whom genotypic resistance testing was routinely available for patients with VF as opposed to a case series or case reports for which the indications for genotypic resistance testing were not reported. Studies that performed routine genotypic resistance testing on all patients with VF provide information on how often PI resistance arises in patients receiving atazanavir. In contrast, the remaining studies were considered likely to be enriched for patients with acquired PI resistance.

Mutations
PI-associated DRMs were defined as those with an HIVDB drug resistance program penalty score for ≥1 PI as of December 31, 2021 [51].

Analyses
The Fisher's Exact Test was used to compare the proportion of each mutation in sequences from patients receiving boosted versus unboosted atazanavir, from patients who were previously ART-naïve versus ART-experienced, and from patients according to whether they had subtype B versus non-subtype B sequences. The Wilcoxon Rank Sum Test was used to compare the median number of mutations between two groups. The Holm's method was used to control for the familywise error rate for multiple hypothesis testing [53]. A binomial regression model was used to examine the relationship between the year of ART initiation and the presence or absence of PI-associated DRMs. To assess the association of covariates with the presence or absence of PI-associated DRMs, a multivariate generalized linear mixed logistic regression analysis was performed using the R package lme4. To account for study heterogeneity, study was included in the model as a random effect.
To identify the patterns of covariation among DRMs and NP-TSMs, we calculated Jaccard similarity coefficients and their standard Z scores for all pair of mutations [54]. To capture conditional dependency among the significantly co-occurring mutation pairs, defined as those pairs that had Jaccard similarity coefficient p < 0.01, we constructed a Bayesian network with a hill-climbing search using the R package bnlearn [55] and created a directed edge network graph using the R package visNetwork [56]. To learn the structure of the Bayesian network of core mutations associated with atazanavir, we excluded sequences containing more than four DRMs in this analysis.
For each sequence containing one or more DRMs, we determined the level of predicted resistance to atazanavir and the levels of predicted cross resistance to lopinavir/r and darunavir/r using the HIVDB drug resistance interpretation system.

Accession Numbers
Sequences in this study had been submitted to GenBank (accession numbers ON058287-ON058987).
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/pathogens11050546/s1, Text S1: The EuResist Network Study Group; Table S1: The complete set of 1497 one-per-person HIV-1 group M protease sequences from persons receiving atazanavir; Table S2: Nonpolymorphic PI treatment-selected mutations (NP-TSMs) occurring in ≥1 sequences from patients receiving boosted or unboosted atazanavir (ATV) as their first PI; Table S3: Studies in PubMed containing sequences from previously PI-naïve patients receiving boosted or unboosted atazanavir (ATV) for which the sequences were not available.