Nationwide Study of Drug Resistance Mutations in HIV-1 Infected Individuals under Antiretroviral Therapy in Brazil

The success of antiretroviral treatment (ART) is threatened by the emergence of drug resistance mutations (DRM). Since Brazil presents the largest number of people living with HIV (PLWH) in South America we aimed at understanding the dynamics of DRM in this country. We analyzed a total of 20,226 HIV-1 sequences collected from PLWH undergoing ART between 2008–2017. Results show a mild decline of DRM over the years but an increase of the K65R reverse transcriptase mutation from 2.23% to 12.11%. This increase gradually occurred following alterations in the ART regimens replacing zidovudine (AZT) with tenofovir (TDF). PLWH harboring the K65R had significantly higher viral loads than those without this mutation (p < 0.001). Among the two most prevalent HIV-1 subtypes (B and C) there was a significant (p < 0.001) association of K65R with subtype C (11.26%) when compared with subtype B (9.27%). Nonetheless, evidence for K65R transmission in Brazil was found both for C and B subtypes. Additionally, artificial neural network-based immunoinformatic predictions suggest that K65R could enhance viral recognition by HLA-B27 that has relatively low prevalence in the Brazilian population. Overall, the results suggest that tenofovir-based regimens need to be carefully monitored particularly in settings with subtype C and specific HLA profiles.


Introduction
According to the Joint United Nations Programme on HIV/AIDS (UNAIDS), 690,000 individuals died in 2019 of human immunodeficiency virus (HIV) infection-related causes. Although the global number of new infections (1.7 million) has been reduced by 23%, since 2010, Latin America has presented a 21% increase. Particularly, Brazil experienced an increase of 17%, with 48,000 new reported infections and 14,000 acquired immunodeficiency syndrome (AIDS)-related deaths, in the same year [1].
Antiretroviral therapy (ART) significantly improves patients' survival, increasing life expectancy and quality [2][3][4]. Furthermore, it plays a relevant role in the prevention of HIV-1 transmission in the population [5][6][7][8]. However, the emergence of drug resistance mutations (DRM) represents a major threat for the continued control of HIV replication, and consequent potential increase in the transmission of viral strains with DRM [9,10]. The prevalence of DRM can greatly vary between different geographical areas, depending

Prevalence of Surveillance Drug-Resistance Mutations
We evaluated the presence of surveillance drug-resistance mutations (SDRM) study population, including SDRM found isolated or in combinations in the same and uncovered an overall prevalence of 84.10% (n = 17,011). SDRM were mostly fou the RT (83.24%; n = 16,836) and included 13,845 sequences with NRTI DRM and sequences with NNRTI DRM. The most common SDRM (Figure 1

A Gradual Alteration on the Prevalence of Surveillance Drug-Resistance Mutations
To investigate the temporal dynamics and evolution of the SDRM landscape in the study population we evaluated yearly SDRM prevalence from 2008 until 2017. The most common mutation across the years was M184V (Figure 1), notwithstanding a decreasing trend and a significant decrease between 2015 (68.29%, n = 1874) and 2016 (53.21%, n = 2684) (p < 0.001) ( Table 2). Similarly, all the other 10 most frequent SDRM followed a decreasing trend along the years with the remarkable exception for K65R and K103N (Figure 1). K103N remained stable in the studied period. While the prevalence of viruses harboring K103N increased between 2014 (40.97%, n = 1180) and 2015 (47.12%, n = 1293; p < 0.001) this was accompanied with a significant decrease in 2016 (42.53%, n = 2145; p < 0.001) and no significant difference between 2008 and 2017 (p = 0.919). Importantly, K65R showed a clear rise along the years (Figure 1). In 2008, the mutation was relatively rare and only found in 2.23% (n = 8) of the population. Subsequently, this mutation became more prevalent and increased significantly between 2013 and 2015 (p < 0.001). Although its prevalence suffered a decrease between 2015 and 2016 (p = 0.0124), K65R prevailed the third most common SDRM in 2016 and 2017 ( Figure 1). Overall, the data strongly supports that, in contrast with other prevalent SDRM, K65R prevalence was increasing in Brazil, raising additional concerns about the efficacy of some of the ART combinations used. In line with this hypothesis was the finding that viral loads for PLWH with K65R were significantly higher than for the rest of the cohort (180,739.60 ± 497,000.10 cop/mL vs. 95,358.26 ± 338,930.67 cop/mL, respectively; p < 0.001). Table 2. Yearly distribution of the most prevalent surveillance drug-resistance mutations in HIV-1 RT and PR.
Percentage number under parentheses. * A statistical difference was found with data from the previous year (p < 0.05). ** Data for 2017 was only available Jan-Apr.

A Shift on Treatment Scheme during the Years
Having observed the clear increase in K65R along the studied yrs, we then decided to evaluate the ART combination schemes in use during the studied period ( Figure 2), aiming at understanding ART usage impact on the SDRM. In 2008, from a total of 359 patients, the most prevalent ART combination was 3TC/AZT/EFV (20.06%, n = 72), followed by 3TC/AZT/LPV (8.36%, n = 30), prevailing the most common treatment till 2014 (25.03%, n = 721) among 2880 patients. From 2015 forward, 3TC/EFV/TDF scheme became the most prevalent. In this year (n = 2744), 24.85% (682) individuals were under 3TC/EFV/TDF drug combination and 19.97% (n = 548) under 3TC/AZT/EFV. In 2017, 3TC/EFV/TDF combination was used on 40.10% (n = 543) of the individuals on ART (n = 1354).

Evidence for Transmission of K65R
To investigate if the transmission of virus harboring K65R could have contributed to the increase in the prevalence of this mutation in Brazil, we performed a phylogenetic analysis on a subset of viral sequences including all the sequences with the K65R mutation and closely related sequences from public databases ( Supplementary Figures S2 and S3). The definition of transmission clusters was performed as previously [37] considering tree branch statistical support and mean genetic distance criteria. We found K65R in at least two HIV-1 sequences in 21 well-delimited transmission clusters (16 from subtype B and 5 from subtype C) (Supplementary Table S1). The size of the inferred transmission clusters ranged from three to 16 individuals. The minimum spanning network analysis of these clusters (Figure 4) supports the occurrence of events of K65R transmission in the study cohort. Furthermore, the analysis of the clinical records from the individuals likely involved in K65R transmission supports this possibility for several cases showing a coincidence of geographic location, proximity in diagnostic dates and similar reported transmission routes (Supplementary Table S1). Despite the significant increased prevalence of the K65R in subtype C (11.26%, n = 347) viruses when compared with B (9.27%, n= 1305) subtype (p < 0.001, Supplementary Figure S2) we found probable events of K65R transmission also in subtype B viruses.
The definition of transmission clusters was performed as previously [37] considering tree branch statistical support and mean genetic distance criteria. We found K65R in at least two HIV-1 sequences in 21 well-delimited transmission clusters (16 from subtype B and 5 from subtype C) (Supplementary Table S1). The size of the inferred transmission clusters ranged from three to 16 individuals. The minimum spanning network analysis of these clusters (Figure 4) supports the occurrence of events of K65R transmission in the study cohort. Furthermore, the analysis of the clinical records from the individuals likely involved in K65R transmission supports this possibility for several cases showing a coincidence of geographic location, proximity in diagnostic dates and similar reported transmission routes (Supplementary Table S1). Despite the significant increased prevalence of the K65R in subtype C (11.26%, n = 347) viruses when compared with B (9.27%, n= 1305) subtype (p < 0.001, Supplementary Figure S2) we found probable events of K65R transmission also in subtype B viruses. We then used an immunoinformatic approach to investigate the predicted impact of K65R on immune-driven selective pressures. We performed artificial neural networkbased predictions of binding to the globally most frequent class I HLAs of all the peptide sequences overlapping K65R and ranging from eight to 12 amino acids with the mutant or the wild-type residue. Interestingly, we found that K65R overlapping peptides were predicted to be recognized by two HLAs (HLA-A03, HLA-B58) in the wildtype version and by one additional HLA (HLA-B27) in the mutant version (Table 3). HLA-B27 has the lowest prevalence (2.23%) in the Brazilian population of all the HLAs tested. Overall, our results suggest that the presence of K65R mutation might contribute for the immune control of the virus in individuals harboring HLA-B27 possibly decreasing its transmission in populations with high prevalence of this HLA. We then used an immunoinformatic approach to investigate the predicted impact of K65R on immune-driven selective pressures. We performed artificial neural networkbased predictions of binding to the globally most frequent class I HLAs of all the peptide sequences overlapping K65R and ranging from eight to 12 amino acids with the mutant or the wild-type residue. Interestingly, we found that K65R overlapping peptides were predicted to be recognized by two HLAs (HLA-A03, HLA-B58) in the wildtype version and by one additional HLA (HLA-B27) in the mutant version (Table 3). HLA-B27 has the lowest prevalence (2.23%) in the Brazilian population of all the HLAs tested. Overall, our results suggest that the presence of K65R mutation might contribute for the immune control of the virus in individuals harboring HLA-B27 possibly decreasing its transmission in populations with high prevalence of this HLA.

Discussion
Acquired and transmitted DRM remain one of the major obstacles towards ART efficacy and AIDS treatment. In case of Brazil, previous regional studies showed distinct moderate patters of transmitted DRM prevalence [11,16,20,27,38], and high DRM rates in patients receiving ART [17][18][19]27,39]. In this study, we aimed at understanding SDRM [26] prevalence across the 27 Brazilian federative units focusing on 20,226 HIV-1 infected patients under different ART schemes, with viral sequence genotyped between 2008 and 2017. In accordance to the Brazilian HIV-1 epidemiology report [36], our cohort was composed mainly by male individuals (55.7%) and the majority of the patients' age varied between 30 and 49 years old (61.44%).
A high prevalence of SDRM was observed with 84.1% of the infected individuals presenting at least one drug resistance mutation. Nonetheless, 68.5% of the genotyped sequences had a SDRM to NRTIs, a value that represents a decrease to what was previously reported in the literature [17][18][19]39,40], following a trend that has been described by Duani et al. and Diaz et al. in previous years [17,40]. However, NNRTI resistance rates (57.95%) did not follow the increasing trend reported by the same authors, presenting a similar value to what has been found in previous works [18,40]. PI resistance mutations (24.82%) also seem to follow the reported decreasing trend [17,40].
Also in line with previous studies is the finding of M184V (65.53%), K103N (40.20%), and M41L (17.21%) as the most common SDRM [17,19,39,40]. M184V, the most prevalent, is a NRTI mutation selected by 3TC and associated with impaired viral fitness and hypersensitization to other NRTI, including AZT and TDF [41][42][43]. Along with the most common mutations, M184V rates presented a decrease over the yrs (76.88% in 2008 to 51.48% in 2017). This decrease on acquired DRM was also observed in other countries and was hypothesized to be due to an improvement on treatment efficacy and also an increased accessibility to ART which enhances patient adherence to therapy [44][45][46]. Moreover, in Brazil, pretreatment genotyping has been extended to children living with HIV under 12 years of age, in case of Mtb co-infection, to pregnant woman and new diagnoses with a sexual partner on ART, since 2013 [47], which may had an impact the on the SDRM general decrease trend.
Nonetheless, it is important to highlight that the percentage of PLWH infected with a virus harboring M148V remained very high and above 50% after 2013. Moreover, K103N and P225H, both frequent NNRTI SDRM, remained stable over the yrs. K103N was previously reported as one of the most commonly acquired SDRM in Brazil [17,18,39,40] in association with the use of EFV [48]. P225H is also selected by EFV and generally occurs in the presence of K103N, synergically increasing EFV resistance [48][49][50]. As for K65R, it is selected by TDF and by other NRTI including abacavir (ABC) and 3TC, decreasing viral susceptibility to most of these drugs [29,31,51].
Previous work reported that K65R mutation presented varied prevalence distribution throughout time in different Brazilian regions [17,19,27]. Strikingly, in our nationwide work in Brazil we found that K65R was the only SDRM that followed an accentuated increase in prevalence over the studied years (2.23% in 2008 to 12.11% in 2017). The increasing and preferential usage of TDF in the clinical practice [21], including in a context of a failing regimen, could be the primordial reason for the significant expansion of K65R as other studies show a higher prevalence of this mutation in patients failing ART treatment [15,[52][53][54][55]. However, although an increasing trend was identified in a cohort from India [56], several studies report a decrease on K65R levels over the years in other countries [48,[57][58][59][60]. Reinheimer and colleagues suggest that the decline on K65R prevalence rates is linked with the increasing usage of single tablet ART regimens and the inclusion of AZT on the used treatment regimen, which has been liked with K65R development suppression [61]. Moreover, Theys et al. work [59] suggests that the found decline of tenofovir K65R selection rate in Portugal, between 2002 and 2010, is mainly caused by changes on treatment guidelines over the years and by the increased usage of combination of TDF and emtricitabine (FTC) [59].
Collectively, these studies suggest that the genetic or sociodemographic characteristics of the population treated with TDF might be influencing the K65R levels.
In 2019, 15.074 individuals were reported to use PrEP in Brazil [1]. Although PrEP consists of a co-formulation of FTC and TDF [62], evidence suggests that the risk of selection for TDF resistance is low and is more likely to occur in cases of undiagnosed HIV infection [63,64]. However, in these cases, high prevalence rates of M184V, which is also selected by FTC, and K65R as found in our study population, might compromise the efficacy of PrEP [65,66].
K65R has been associated with diminished viral fitness and replication, which has been linked with decreased transmission capacity [67,68]. Even though several studies support almost inexistent K65R transmission rates [11,69,70], Rhee and co-works present evidence of higher levels of transmitted K65R, especially in low-and middle-income countries [71]. In our phylogenetic analysis we found 21 transmission clusters containing at least two sequences with the K65R mutation. Moreover, minimum spanning network analysis supported the occurrence of events of K65R transmission in at least seven of these clusters. This might be suggestive that the efficacy of PrEP might be compromised in this setting. Indeed, we found much higher viral loads in patients harboring K65R mutation when compared with the rest of our cohort population. Similar viral loads were observed between individuals with and without tenofovir resistance by The TenoRes Study Group [72], which can be associated with increased mutation transmission. Furthermore, our results suggest that K65R overlapping peptides can be increasingly recognized by HLA-B27, a HLA that is linked with lower viral loads and slower diseased progression to AIDS [73,74]. The fact that Brazil presents a low prevalence of the HLA-B27, when in comparison with other regions [75][76][77], might have contributed to K65R expansion in this country. However, it is relevant to point-out that our results are based on in-silico prediction and not taking in consideration the possible impact on immune responses of the combination of different mutations in the same virus. Thus, it is relevant to address this topic in the future by performing functional studies.
Despite of the observed decreasing trend in the prevalence of some major SDRM in the studied Brazilian cohort, a high prevalence of these mutations was still verified. Indeed, our results show that the alteration performed in 2013 enlarging the criteria to include more PLWH in the recommendation for baseline SDRM testing had a positive impact but was still insufficient. We propose that the lack of universal baseline HIV-1 DRM screening to inform on effective ART regimens resulted in high levels of SDRM, such as M184V, K103N, and M41L underlying many cases of treatment failure in Brazil not only from 2008-2012 but also continuing from 2013-2017. Furthermore, we observed a clear increase in K65R reverse transcriptase SDRM that is an additional problem in Brazil that could have been aggravated by the circulation of HIV-1 subtype C and the HLA class I makeup of the population. This increase in K65R was mainly found in association with the use of TDF and is particularly relevant in combination with the high M184V levels found in the study population suggesting that the efficacy of PrEP might be compromised.
Overall, our results support that some of the drugs most frequently used in Brazil might be compromised due to the high frequency of SDRM and that baseline drug resistance testing should be universal and mandatory as it is the best way to promote personalized selection of the most optimized ART regiment.

Study Population
The collection of the patient's data was performed anonymously after approval by the Brazilian national ethic committee through the protocol CAAE 53757016.0.0000.5504. All records from HIV-1-infected individuals followed by the Specialized Assistance Services on Sexually Transmissible Diseases and HIV/AIDS obtained at all 27 Brazilian federative units from 01/01/2008 to 04/30/2017 were curated in a relational database. A total of 20,226 HIV-1 infected patients were selected for this study according to the following inclusion criteria: (i) be under ART treatment; (ii) have at least one associated partial HIV-1 genome sequence available (complete HIV protease (n = 20,223) and complete reverse transcriptase (n = 20,214)). HIV-1 sequencing was performed as part of the routine clinical testing with commercially available HIV-1 genotyping systems based on Sanger sequencing.

HIV-1 Subtypes and Drug Resistance Mutations
HIV-1 subtypes were identified by relying on the consensus from least two of the three utilized subtyping tools: SNAPPY [78], jpHMM [79] and Stanford HIV Drug Resistance database (https://hivdb.stanford.edu/, accessed on 1 May 2021). All sequences were also analyzed using the Stanford HIVdb Genotypic Resistance Interpretation Algorithm to evaluate and interpret the presence of DRM in each sequence.

Phylogenetic Analysis of K65R Sequences and Transmission Clusters
All sequences having the K65R were separated by subtype and used to query local and public databases to identify highly related HIV-1 sequences. The resulting sequences from the most common subtypes B and C were aligned with MAFFT v7.309 [80] and used to make a phylogenetic reconstruction using PhyML v3.0 [37]. The best fitting substitution model was GTR + G4 + I, determined by PhyML SMS using AIC [81]. The heuristic trees search was performed using SPR methods. Bayesian evolutionary analyses were performed using BEAST v1.10.4 [82,83] with GTR + G4 + I, as the nucleotide substitution model. This phylogenetic representation was used to infer the transmission clusters as previously [84]. Minimum spanning network analysis of the sequences identified in transmission clusters with more than one K65R mutant was performed with PHYLOViZ [85].

HLA Prevalence in the Brazilian Population
The Allele Frequencies Database [87] was used as the source of this information for the Brazilian population. To obtain cohesive and yet representative nationwide results, the bone marrow registry (REDOME) data was selected. Since this data is separated by Brazilian state, we calculated the proportional values for a national level estimation.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ijms22105304/s1, Figure S1: Presence of SDMR among subtypes, Figure S2: Phylogenetic representation of a subset of HIV-1 subtype C sequences isolated in Brazil with the K65R mutation and closely related sequences from databases, Figure S3: Phylogenetic representation of a subset of HIV-1 subtype B sequences isolated in Brazil with the K65R mutation and closely related sequences from databases. Table S1: Phylogeny inferred HIV-1 transmission clusters with more than one sequences harboring the K65R mutation, Table S2 Funding: This work has been funded by Portuguese National funds, through the Foundation for Science and Technology (FCT) (project UIDB/50026/2020 and UIDP/50026/2020; by the projects NORTE-01-0145-FEDER-000013 and NORTE-01-0145-FEDER-000023, supported by Norte Portugal Regional Operational Programme (NORTE 2020), under the PORTUGAL 2020 Partnership Agreement, through the European Regional Development Fund (ERDF) and by Gilead Génese PGG/009/2017. ASP and PMMA were funded by FCT PhD scholarships PD/BD/127827/2016, and PDE/BDE/113599/2015, respectively.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Brazilian National Ethics Committee through the protocol CAAE 53757016.0.0000.5504.
Informed Consent Statement: Patient consent was waived since all data was anonymous and collected retrospectively.

Data Availability Statement:
The HIV-1 sequences selected for this study were made available in GenBank (accession numbers will be included here).