Molecular evolution of HIV-1 integrase during the 20 years prior to the first approval of integrase inhibitors

Detailed knowledge of the evolutionary potential of polymorphic sites in a viral protein is important for understanding the development of drug resistance in the presence of an inhibitor. We therefore set out to analyse the molecular evolution of the HIV-1 subtype B integrase at the inter-patient level in Germany during a 20-year period prior to the first introduction of integrase strand inhibitors (INSTIs). We determined 337 HIV-1 integrase subtype B sequences (amino acids 1–278) from stored plasma samples of antiretroviral treatment-naïve individuals newly diagnosed with HIV-1 between 1986 and 2006. Shannon entropy was calculated to determine the variability at each amino acid position. Time trends in the frequency of amino acid variants were identified by linear regression. Direct coupling analysis was applied to detect covarying sites. Twenty-two time trends in the frequency of amino acid variants demonstrated either single amino acid exchanges or variation in the degree of polymorphy. Covariation was observed for 17 amino acid variants with a temporal trend. Some minor INSTI resistance mutations (T124A, V151I, K156 N, T206S, S230 N) and some INSTI-selected mutations (M50I, L101I, T122I, T124 N, T125A, M154I, G193E, V201I) were identified at overall frequencies >5%. Among these, the frequencies of L101I, T122I, and V201I increased over time, whereas the frequency of M154I decreased. Moreover, L101I, T122I, T124A, T125A, M154I, and V201I covaried with non-resistance-associated variants. Time-trending, covarying polymorphisms indicate that long-term evolutionary changes of the HIV-1 integrase involve defined clusters of possibly structurally or functionally associated sites independent of selective pressure through INSTIs at the inter-patient level. Linkage between polymorphic resistance- and non-resistance-associated sites can impact the selection of INSTI resistance mutations in complex ways. Identification of these sites can help in improving genotypic resistance assays, resistance prediction algorithms, and the development of new integrase inhibitors.


Background
The HIV-1 integrase catalyses the integration of the reverse transcribed viral DNA into the host genomic DNA via a two-step process. In its active form the integrase forms a tetramer. The monomeric enzyme consists of 288 amino acids (aa) and contains three functional domains: the N-terminal zinc-binding domain (NTD, aa , the central catalytic core domain (CCD, aa , and the C-terminal DNAbinding domain (CTD, aa 213-288) [1][2][3][4]. Each region comprises motifs essential for the proper function of the enzyme, e.g. the zinc finger motif H12-H16-C40-C43 in the NTD, the active site D64-D116-E153 in the CCD, and the minimal nonspecific DNA-binding region ranging from I220 to D270 in the CTD [2][3][4][5][6].
Raltegravir was the first integrase strand inhibitor (INSTI) to be approved in Europe in 2007, followed by elvitegravir in 2012 and dolutegravir in 2014. Bictegravir [7] and cabotegravir [8] are in clinical trial development.
INSTIs target the CCD, thereby inhibiting the strand transfer of the double-stranded viral DNA into the host genome [1]. Various allosteric inhibitors of integrase (ALLINIs), which modulate integrase multimerisation [9] and interfere with the cellular transcription factor LEDGF/p75 [10] are in development but have not so far made it further than Phase I clinical trial [11].
Due to its high rate of replication, mutation, and recombination, HIV is a virus of high genetic variability. The viability of virus variants in turn is limited by structural and functional constraints. At the same time, variants in the viral quasispecies can be selected by the pressure of the human immune system [12,13] or antiretroviral treatment (ART) [14,15]. In general, the HIV diversity at the intra-patient level increases during the course of infection [16] driven by both drift and selection [17]. During transmission to a new host, several stochastic and selective bottlenecks reduce the viral diversity to a few variants [18], and factors that contribute to shaping the HIV diversity at the inter-patient level are extensively discussed [19][20][21][22][23][24].
Naturally occurring polymorphisms can affect the genetic barrier to drug resistance by influencing the selection of resistance mutations, enzymatic activity, and replicative capacity [22,25,26]. Epistatic interactions between polymorphisms can further modulate viral fitness and the development of drug resistance in complex ways and have been shown to play an important role in the HIV-1 protease and reverse transcriptase [27][28][29][30]. Thus, to understand the selection of resistant variants in the presence of INSTIs, it is important to investigate the evolutionary dynamics of the polymorphic sites in the integrase.
The prevalence of HIV-1 integrase polymorphisms and INSTI resistance mutations has been investigated before in INSTI-naïve individuals [26,[31][32][33][34][35] and ART-naïve individuals [35][36][37][38][39][40]. However, time trends and covariation of complex mutation patterns preceding the availability of INSTIs have not so far been analysed. The aim of this study was to investigate covarying clusters of naturally occurring resistance-and non-resistanceassociated amino acid variants and their frequencies over time at the inter-patient level to consider their potential relevance for INSTI resistance. To this end, HIV-1 integrase sequences were obtained from samples of ART-naïve individuals newly diagnosed with HIV-1 between 1986 and 2006, a 20-year period prior to the first approval of INSTI in Germany [41].

Study population
Plasma samples from individuals newly diagnosed with HIV-1 between 1986 and 1996 (N = 167) were archived at the former diagnostic unit of the National AIDS Centre in Germany and stored at −40°C. The date of HIV-1 diagnosis is the same as the plasma sampling date for the HIV-1 integrase genotyping. Plasma samples from individuals newly diagnosed with HIV-1 between 1997 and 2006 (N = 170) were collected for the German HIV-1 Seroconverter Study [42][43][44][45] and stored at −70°C . These plasma samples were taken within 12 months after diagnosis. In total, the study population comprised 337 individuals (Table 1).

Assessment of phylogenetic bias
A unique anonymising code provided with the mandatory report of newly HIV-1 infected cases to the Robert Koch-Institute ensured that only one HIV-1 sequence per patient was included in our dataset. However, to investigate whether our results were biased by the overrepresentation of phylogenetically closely related sequences (i.e. sequences originating from direct transmission events) we identified clusters of sequences with a very small phylogenetic distance and replaced them with one representative sequence (the sequence from the patient who was diagnosed first). For clustering, a multiple sequence alignment (MSA) including an HIV-1 group N reference sequence (Acc. AJ006022) was generated with clustalW (BioEdit, v7.2.5, Tom Hall) and endtrimmed to nt 4230-5064 (Acc. K03455). We computed a maximum likelihood phylogeny with 500 times nonparametric bootstrap with replacement using the program RAxML [49]. The HIV-1 group N sequence was used to root the maximum likelihood tree. Applying the program Transmic [50], a set of phylogenetically closely related sequences was identified if the mean of all pairwise patristic distances did not exceed a threshold of 0.015 expected nucleotide substitutions per site and the most recent common ancestor node had a bootstrap support of 0.9 [50,51]. We detected 11 clusters, each comprising two sequences. The final reduced dataset consisted of 326 sequences. All analyses were then performed with the full dataset and the reduced dataset to assess any potential bias. All reported results were confirmed on the basis of both datasets while all numbers reported in the results section are given for the full dataset.

Amino acid variability
MSAs for the full (337 sequences) and the reduced (326 sequences) dataset were generated with clustalW (BioEdit, v7.2.5, Tom Hall) and end-trimmed to nt 4230-5064 (Acc. K03455). The nucleotide sequences were then translated into amino acid sequences. Nucleotide sequence ambiguities of codons were resolved during translation. An X was assigned if multiple amino acids resulted from the translation of codons containing nucleotide sequence ambiguities to avoid consideration of amino acid variants only present in mixtures due to PCR and sequencing errors. The frequencies of amino acids at positions 1 to 278 of the HIV-1 integrase were calculated. A position was defined to be polymorphic if, on the basis of the full and the reduced dataset at this position, a total of 1% or more amino acid variants were present when compared to the consensus B sequence (http://hivdb.stanford.edu/pages/documentPage/consensus_amino_acid_sequences.html). All other positions were defined as conserved.
We also analysed the Shannon entropy E i , which quantifies the degree of variability at a single amino acid position i [1,278], according to where P(X i = x ) denotes the probability of observing the particular amino acid x at position i [52]. A high value of entropy E i indicates high amino acid variability at position i.

Amino acid time trends
The samples were grouped into five time periods (1986-1989, 1990-1994, 1995-1998, 1999-2002, and 2003-2006) to obtain a more homogenous distribution and to render the time trend analyses more robust (Additional file 1: Table S1). Next, we calculated the frequencies of all amino acids at positions 1 to 278 of the HIV-1 integrase within the periods. We then fitted a linear regression by minimizing the sum of least square deviation between detected and predicted frequencies y and f(a, b) for each period i according to where a denotes the slope and b the intersect of the linear function. The variable x i denotes the distance of the centre of period i [in years] to the centre of the first period. In order to detect significant time trends, we generated 10,000 bootstrap samples by drawing sequences from the original sequence set with replacement. For each resampled set, we computed mutation frequencies and performed the linear regression as described above. Raw P values were then computed from the bootstrap distribution of fitted slopes a in analogy to Katchanov et al. [53], i.e. we computed P þ ¼ #a≤0:001 10000 to test whether the frequency of that mutation is significantly increasing by at least 0.1% per year (H 0 : a ≤ 0.001 vs. H 1 : a > 0.001) and conversely P − ¼ #a≥−0:001 10000 to assess whether it is significantly decreasing by at least 0.1% per year (H 0 : a ≥ − 0.001 vs. H 1 : a < − 0.001). P values (P = min(P − , P + )) were subsequently corrected for multiple testing using the false discovery rate (FDR) method by Benjamini-Hochberg [54]. Time trends were considered significant if they were identified on the basis of the full and the reduced dataset with a FDR corrected P < 0.05.

Amino acid covariation
In order to determine if amino acid positions are covarying, we applied direct coupling analysis where direct correlations are disentangled from transitive correlations. We computed evolutionary coupling terms ec ij using the recently developed plmc tool (https://github.com/debbiemarkslab/plmc) [55]. This tool infers couplings by fitting a Potts model to the MSA using a pseudo likelihood approach with L 2 regularisation. We used default regularisation parameters of λ 1 = 0.01 and λ 2 = 100 for single site contributions and direct couplings terms e ij (α,β), respectively. To correct for a phylogenetic bias by reweighting neighbouring sequences, we chose the default parameter θ = 0.01. Ambiguous amino acids (X) were discarded during inference. As expected, the inferred coupling terms are approximately normally distributed with mean 0 and standard deviation 1/λ 2 .
The direct coupling terms e ij (α,β) describe the direct correlation of two amino acids α and β at positions i and j. The overall interaction of two positions i and j are given by the evolutionary coupling term ec ij , which is the Frobenius norm of the direct coupling terms with average product correction (APC), to suppress phylogenetic bias effects [56]: To identify significant coupling terms, z-scores were computed according to where μ(e * * ) and σ(e * * ) denote the mean and standard deviation of the estimated evolutionary coupling terms between all position pairs i < j. Z-scores were then converted into P values and corrected for multiple testing using the Benjamini-Hochberg method [54]. We used P < 0.005 to detect significant couplings [57], which corresponds to a Bayes factor of 14 typically used in Bayesian inference [58].
To ensure the robustness of the statistical inference, Halsey et al. [59] recently proposed to combine power analysis with P value based statistical analysis. In line with these recommendations, we generated 1000 MSAs by sampling sequences with replacement from the original alignment. For each resampled MSA, we inferred the coupling terms and performed the statistical tests as outlined above. In line with Halsey's recommendations, we only reported coupling terms that were significant in at least 95% of the resamplings, i.e. where the respective test power was >0.95. An interactive plot to explore the detailed results of the direct coupling analysis was generated using the EVZoom tool (https://github.com/debbiemarkslab/ EVzoom) and can be accessed through http://page.mi.fu-berlin.de/msmith/couplings_integrase.html.

Frequency of INSTI resistance mutations
As expected, no major INSTI resistance mutations were detected in this ART-naïve study population from the period prior to INSTI release. Moreover, all major INSTI resistance-associated positions were fully conserved, and amino acids at these positions corresponded to the respective INSTI-sensitive consensus B amino acids. However, some minor INSTI resistance mutations and INSTI-selected mutations were identified (Table 2).
At nine positions (6, 11, 72, 101, 135, 154, 165, 201, 265) an amino acid exchange between two variants was identified. In contrast, at sites 119 and 124 increasing polymorphy was observed since variants decreased in frequency without being replaced by another variant. Finally, variants 122I and 256E increased in frequency without replacing another variant, although a decrease of T122 and D256 was observed that was close to statistical significance (Fig. 3). Variants L101I, T122I, and V201I with increasing frequency as well as M154I with decreasing frequency were associated with INSTI resistance (Fig. 3).

Amino acid covariation
Overall, 42 significant couplings were identified (Table 3). 10 couplings appeared within the NTD, 10 within the CCD, and two within the CTD. 13 couplings were observed between the NTD and the CCD, five between the CCD and the CTD, and two between the NTD and the CTD. In total, 25 and 28 couplings involved at least one position within the NTD and the CCD, respectively (Fig. 4).
17 amino acid variants with temporal trend were found to covary with other time-trending variants. The individual time trends were generally compatible, i.e. both time trends were concordant if the coupling terms were positive and discordant if the coupling terms were negative (Table 4).
We found the highest amino acid variability determined by entropy, time trends, and direct coupling analysis within the CCD, the NTD, and between CCD and NTD. Sites important for enzymatic activity were in general conserved, however, some positions involved in binding the cellular cofactor LEDGF/p75 (sites 125, 165) and within the minimal nonspecific DNA binding region (sites 220, 230, 232, 234, 256, 265) were polymorphic with an overall variability ≥5%. Covariation between positions 125, 165, 256, or 265 and other sites was observed, and the DNAbinding site 234 covaried with the DNA-binding site 253. The most frequent substitutions were T125A, V165I, I220L, S230 N, D232E, L234I, D256E, and A265V. All of them occurred within the same biochemical class of amino acids, with the exception of T125A that represents a switch from a hydrophilic to a hydrophobic amino acid.   [67] as an overlay of protein database structures 1ex4 [3] and 1k6y [4]. The NTD (aa 1-49) is coloured yellow, the CCD (aa 50-212) green, and the CTD (aa 213-288) red. Polymorphic sites are coloured grey The effect of this switch is not known and should be investigated experimentally. A time trend in frequency was observed for variants T125A, D256E, and A265V. The knowledge about the variability of the integrase should be taken into account for the design of genotypic resistance assays.
Most time trends were based on an exchange of two amino acids, however, a general diversification was observed at sites 119 and 124. 17 out of 22 amino acid variants with increasing or decreasing frequency covaried among each other. In general, we observed a concordant time trend for pairs with a positive direct coupling term and a discordant time trend for pairs with a negative direct coupling term. Exceptions to this rule were couplings between positions 154-265 and 201-256 (Table 4). The time trends for the individual variants 154I-A265, M154-265 V, and V201-256E were discordant despite positive direct coupling terms. The reason for this may be unidirectional coupling, i.e. 154I requires the presence of A265, but not vice versa. Likewise, the time trends of M154-A265, 154I-265 V, and 201I-256E were concordant despite negative direct coupling terms.
The prevailing concordance of significant time trends and significant couplings in our study suggests the selection of coevolving epistatic clusters. However, due to the transmission bottlenecks [18], genetic drift may be another viable explanation for the observed time trends in the frequency of certain amino acid variants. The role of genetic drift in HIV evolution is debated and has been quantified to some extent at the level of intra-patient evolution and for known transmission pairs rather than on inter-patient level [60,61]. Genetic drift on interpatient level requires inheritance of selectively neutral substitutions. Large parts of the integrase may be under negative selection to maintain enzymatic functionality; nevertheless, particular positions and certain substitutions of the integrase may be selectively neutral. Therefore, we considered the possibility of genetic drift for all time-trending substitutions by assessing whether there Positions showing trends in INSTI resistance-associated variants are shaded grey. *** time trend is significant at the P < 0.01 level, ** time trend is significant at the P < 0.05 level was evidence for (i) inheritance of the time-trending substitutions and (ii) whether the time-trending substitutions were selectively neutral. Ad (i): We statistically compared the mean patristic distance of random sequences versus the mean patristic distance of sequences carrying a specific time-trending substitution to investigate if time-trending substitutions appeared more frequently in phylogenetically closely related sequences (Additional file 1: Table S1). Sequences carrying the time-trending substitutions 72 V, 154I, 165I, and 265 V had significantly smaller mean patristic distances than random, by this indicating inheritance. Interestingly, all of these substitutions decreased in frequency over time. Ad (ii): First, we performed Tajima's D test [60,61] with a result of D = −1.44, by this indicating negative selection. Next, we calculated the ratio of nonsynonymous over synonymous mutations (dN/dS ratio) [60,61], finding that most regions of the integrase were under strong negative selection, including sites 72, 154, 165, and 265 (Additional file 2: Figure S1). In summary, we could not observe a clear contribution of genetic drift to the time trends of the examined substitutions.
By using Sanger sequencing of bulk RT-PCR products with a sensitivity of approximately 30% [62,63] and excluding ambiguous amino acids in our analyses we could only investigate the major virus variant from each patient sample. Minor variants and linkage between minor variants can only be investigated by using more sensitive techniques like single genome sequencing (SGS) or next generation sequencing (NGS) [63][64][65]. To minimize the probability that technical errors during RT-PCR and Sanger sequencing lead to false positive predictions with regard to coupling terms, we combined our direct coupling analysis with a power analysis, which essentially requires that an amino acid pair has to be present in multiple sequences to be repeatedly identified by direct coupling analysis. Recently, the use of covariation methods as a measure of coevolution has been questioned by Talavera et al. [66]. Based on a computational study, the authors point out that a strong covariation signal is caused by a low evolutionary rate. We therefore assessed our results accordingly but could not find a relation between the rarity of pairwise substitutions and high coupling terms or the occurrence of single substitutions in couplings (Additional file 3: Figure S2).
16 minor INSTI resistance mutations and 11 INSTI-selected mutations were observed as naturally occurring in our ART-naïve study population, which originated from the time prior to INSTI approval. Among these resistance-associated variants, three increased in frequency over time and seven covaried with non-resistance-associated variants. The complex Evolutionary coupling terms ecij between position i and j with P < 0.005 and power > 0.95 are given interdependent evolution of these mutations might control enzymatic activity and replication capacity independent of selective pressure through INSTIs at the inter-patient level. Indeed, accessory drug resistance mutations that compensate viral fitness are often already polymorphic in drug-sensitive HIV-1, suggesting that these mutations may naturally enhance viral fitness and virulence with progression of the HIV-1 epidemic [21,22]. INSTI-independent linkage between non-resistance-associated sites and resistanceassociated sites or sites targeted by INSTIs can affect the selection of resistance mutations in the presence of INSTIs. This knowledge should be taken into account for the improvement of resistance prediction algorithms as well as for the development and preclinical evaluation of new INSTIs and ALLINIs. Deeper analyses of the observed resistance-associated variants are needed to evaluate their clinical relevance. In particular, those with naturally increasing frequencies that were linked to covariation should be investigated, i.e. L101I, T122I, and V201I. The absence of major INSTI resistance mutations in our ART-naïve study population underscores the suitability of INSTIs for first-line antiretroviral regimens.
Because the analysed dataset was rather small (n = 337), our results may require further validation from the analysis of larger, independent datasets. Due to the relatively small number of samples, some of our results might not have reached statistical significance, e.g. the temporal trend of T122 and D256. Generally, given the small sample size, overrepresentation of almost identical sequences (i.e. from transmission chains) may profoundly bias any downstream analysis of time trends and covariation patterns. To assess whether our analyses were affected by such sampling bias, we additionally performed them using a reduced dataset in which clusters of closely related sequences were replaced by one representative only. The results obtained from the reduced dataset confirmed all results obtained from the full dataset.

Conclusions
Our aim was to analyse the molecular evolution of the HIV-1 integrase prior to the approval of INSTIs and, thus, INSTI selective pressure at the inter-patient level. We found significant time trends in the frequency of certain amino acid variants, suggesting ongoing adaptation of the enzyme. Upon closer inspection, we found that amino acid variants with significant time trends covaried with other time-trending variants. Such a linkage may impose constraints that determine the evolutionary

Additional files
Additional file 1: Table S1. Statistics assessing whether time-trending substitution could have resulted from inheritance. (DOCX 11 kb) Additional file 2: Figure S1.