HIV-1 replicates via a low-fidelity polymerase with a high mutation rate; strong conservation of individual nucleotides is highly indicative of the presence of critical structural or functional properties. Identifying such conservation can reveal novel insights into viral behaviour. We analysed 3651 publicly available sequences for the presence of nucleic acid conservation beyond that required by amino acid constraints, using a novel scale-free method that identifies regions of outlying score together with a codon scoring algorithm. Sequences with outlying score were further analysed using an algorithm for producing local RNA folds whilst accounting for alignment properties. 11 different conserved regions were identified, some corresponding to well-known cis-acting functions of the HIV-1 genome but also others whose conservation has not previously been noted. We identify rational causes for many of these, including cis functions, possible additional reading frame usage, a plausible mechanism by which the central polypurine tract primes second-strand DNA synthesis and a conformational stabilising function of a region at the 5′ end of env.
Human Immunodeficiency Virus (HIV) infection remains a significant global health burden, with an estimated 36. 7 million people worldwide living with the virus and 1. 0 million AIDS-related deaths in 2016 [1]. Despite the success of existing interventions in reducing the annual incidence of HIV, there is still substantial scope for further progress in limiting new infections and optimising diagnosis and treatment of those with HIV. Improving our knowledge of the viral genomic structure and function permits better understanding of the viral lifecycle and host/virus interactions and may suggest interventions to target viral replication and survival. The HIV-1 genome is one of the most intensively studied genetic sequences. It contains a large number of cis-acting regions whose function depends on the structures into which the RNA folds. Many of these have been studied and solved at a secondary structure level and for some there are three-dimensional data. The whole genome has been analysed biochemically at a secondary structure level and many of the known functional regions have been mapped [2]. Much of the genome, however, has been determined to possess conserved structure by such techniques, but with no functional data. Similarly, there are regions where the function of the individual sequence nucleotides, either directly or through coding function, is the main constraint on the sequence and primary structure dominates over secondary structural requirements. HIV-1 has a high mutation rate, with genome-wide estimates varying from 5. 3 × 10−5 to 5. 8 × 10−3 per base per cycle, values of estimates depending upon a number of factors including how reverse transcriptase variability is measured, estimates of recombination rates, and the contribution of cellular apolipoprotein B messenger RNA-editing enzyme-catalytic polypeptide-like 3 (APOBEC3) proteins to variability [3–7]. Notwithstanding variations in the mutation rate from base to base, this number may be sufficiently high to expect every possible point mutation to occur every day in an individual infected with HIV-1 [8]. Conservation within observed sequences is therefore significant and strongly indicative of evolutionary selection pressure, and points towards the presence of structures or functional regions necessary for viral survival. Detection of such conservation can help direct experimental searches to reveal novel structural and/or functional aspects. Without a priori knowledge of the dimensions of individual structural or functional elements, analytical techniques need to be able to detect conserved elements regardless of scale. To investigate further the possible structural form of the viral RNA and to seek evidence of conservation within the genome that might reveal previously unidentified cis-acting regions and functions, we performed an intensive bioinformatic analysis of HIV-1 genetic data. We analysed sequence variation, controlling for codon usage, using a novel scale-free method we have developed for finding regions of conservation when the underlying driver of conservation and hence the order of magnitude of the size of conserved regions is not known. This method differs from previous analyses in that it makes no assumptions about the scale of features of interest. As such it will reveal regions that are missed by previous methodologies. Our analysis reveals eleven regions of sequence conservation. Some of these coincide with clusters of established important cis-acting functions such as splice sites and their regulatory regions but for others there is no existing explanation; however, some highly likely causes can be attributed to a number of these. Analysing the secondary structure of these sequences reveals some known, but also a number of novel highly conserved structures including a previously undescribed highly stable structure around the central polypurine tract, with clear parallels between this and the structure we identify for the 3′ polypurine tract. Sequence preserved regions such as these represent promising targets for design of HIV therapeutics.
For each gene in the HIV-1 genome, as well as for the entire genome, non-recombinant B-subtype sequence data collated between 2009 and 2015 were downloaded from the Los Alamos National Laboratory Filtered Web Alignments database [9]. This database ensures only one sequence is included per patient, that a single representative is included of very similar sequences, and that sequences unlikely to represent natural viable viruses are excluded. Sequence data were further filtered to exclude sequences with uncertain or missing nucleotides in the region being studied, and to exclude sequences whose length indicated that they were incomplete or contained large insertions. A full listing of the sequences used can be found in the supporting information (S1–S10 Tables). Sequences for each gene were aligned at amino acid level using MUSCLE [10]. Each gene was analysed at codon level with Mathematica [11] using the normalised mean pairwise distance (nMPD) method, as previously described [12], subject to the modification that positions with invariant tryptophan and methionine codons were recorded as uninformative rather than having nMPD set to 1. The nMPD method scores each codon locus by summing the pairwise Hamming distances between all codons. It then normalises to take into account amino acid usage and codon bias, by dividing by the expected distance sum taking into account the amino acids used at that locus and the distribution of codons encoding those amino acids across the genome. At each codon locus a non-negative nMPD score is produced, with low scores representing unexpected conservation of nucleotides after amino acid conservation and codon bias have been taken into account. Such unexpected conservation may represent important structural elements or important functional elements (or both). The analysis of HIV alignments using the nMPD method raised an additional issue that the gap percentage is higher than that found in alignments of pathogens previously studied using this method, meaning that the way in which gaps are handled can have a substantial effect upon nMPD scores and consequently on the presence or absence of runs of low scores. Gaps were handled by defining the Hamming distance between a gap and any other nucleotide (including a gap) to be zero, but by labelling as uninformative codons where the gap percentage level exceeded 10%. This choice ensures that a minority variant in which a small motif is absent from or present in an otherwise highly conserved region does not result in the entire region being excluded from recognition. The original nMPD analysis [12] used sliding windows to identify runs of unexpectedly low scores, corresponding to regions of amino acid conservation. Here we instead used the novel scale-free analysis method we recently described [13]. This method is more versatile in that it more easily allows agnostic treatment of uninformative regions and it makes no a priori assumptions about the length of conserved motifs, although it still requires conserved regions to be relatively contiguous. The method converts nMPD data to have the appearance of a series of random walks, where at each stage of analysis the presence of an unusually steep descent corresponds to a region of unexpected nucleotide conservation. Regions determined by the above algorithm to be significantly conserved were aligned to the HXB2 ([14], GenBank accession K03455) and NL4-3 ([15], GenBank accession AF324493) sequences for reporting and analysis. Regions determined by the algorithm to be significantly conserved were further analysed by taking the relevant section of the whole-genome alignment and folding the region using RNAalifold [16,17], with default parameters except for using ribosum scoring and disallowing helices of length 1 (lonely pairs). RNAalifold combines phylogenetic information from an input sequence alignment with thermodynamic parameter calculation to determine a consensus secondary structure. The output allows local structures to be analysed for conserved features of interest and conserved loci and reciprocal mutations. Structures are presented in their original RNAalifold format and not modified on the basis of other published structural data as both our novel nMPD analysis and RNAalifold are designed to generate data from multiple sequences, whereas structural prediction from biochemical analyses such as selective 2′-hydroxyl acylation analyzed by primer extension (SHAPE) use only a single sequence as input. Generating a hybrid approach to constrain secondary structure using primary structure data from multiple sequences and e. g. SHAPE data from a single separate sequence would be of uncertain validity. As an additional validation step, we acquired sequences and undertook the nMPD analysis step for A-subtype sequence data (including sequences labelled with sub-subtypes within the A-subtype), in the same way as described for the B-subtype sequence data, mutatis mutandis. A full listing of the sequences used can be found in the supporting information (S11–S20 Tables). The purpose of this step was to provide additional validation of the main B-subtype nMPD analysis results and we did not perform further structural analysis of the A-subtype data.
nMPD analysis was performed for each gene, using between 309 and 1236 sequences (in total 3651 distinct sequences were used). The output of nMPD analysis is shown in S1–S9 Figs. Five genes (pol, tat, vpu, env, nef) contained regions that were deemed significantly conserved (Fig 1 and Table 1 for summaries of regions and features of interest found; Figs 2–4 for plots of cumulative analysis for genes in which significantly conserved regions were found; plots of cumulative analysis for remaining genes are in S10–S13 Figs). The conserved region found first in env, and the fourth region identified reading from its 5′ end—HXB2 nucleotide reference 7662–8051 (NL4-3 7652–8041) —corresponds well with the well-known highly conserved Rev-response element (RRE) [21–23]. Our method clearly picks up how strongly conserved this region is and delineates its location (see Table 1 for a summary of the statistics, and Fig 2 for a visual representation). Successful detection of this region serves as a positive control for our method. A local fold of this region is shown in Fig 5. The RNAalifold representation is very similar, albeit not identical, to other predictions for folding of the RRE generated by other methods (e. g. reference [2]). The conserved region found by the first iteration of our algorithm in pol, and the second region identified reading from its 5′ end (Fig 2) —HXB2 nucleotide reference 4749–5093 (NL4-3 4749–5093) —contains the central polypurine tract and associated G-rich sequences suggested to be important for RNA dimer stabilisation via G-quartet formation [24]. The region also contains a set of structured stem-loops that contain the A1 splice acceptor site and D2 splice donor site, plus associated exon splicing enhancer elements: ESEVif [25], ESEM [26], and a novel exon splicing enhancer surrounded by the ESEM1 and ESEM2 domains, which was initially identified with the HEXplorer algorithm [27]. The region also contains a splice silencing element, G4 [25], which in our fold of the region is located in the stem of the structure containing the D2 splice site (Fig 6), and a postulated regulatory G run, GI2-1 [28], located before the vif start site. There is an enhancer motif present between the G4 element and the silencing regulatory G run, termed ESE2b, and directly upstream of this the recently described ESS2b domain; mutations in or changes in accessibility of these domains have been shown to change the proportions of spliced products and affect viral infectivity [29–32]. The region finishes just 3′ of the rarely used D2b splice donor site. In the region is a motif previously identified as evolutionarily conserved, termed stem-loop containing splice acceptor 1 (SLSA1) [33,34], and although we note that the computational structural prediction by the method we have used does not reproduce the stem-loop of the original prediction, this does not detract from the presence of the conserved nucleic acid motif. The region also contains the start of vif (i. e. there are overlapping open reading frames, resulting in conservation of the wobble position nucleotides), and the region ends just 5′ of the pol stop codon. A local fold of the region is shown in Fig 6. One region of interest is highlighted in vpu (Fig 4), HXB2 nucleotide reference 6206–6226 (NL4-3 6202–6222), just 5′ to the start of env. The first region of interest highlighted in env (Fig 2) is nearby, at HXB2 nucleotide reference 6324–6611 (NL4-3 6320–6607). Not far 5′ of these regions is the region of interest highlighted in tat (Fig 3), HXB2 nucleotide reference 5951–5998 (NL4-3 5950–5997), containing the A4a, A4b, A5, A5a and A5b splice sites, in which possible conserved stem-loop structures have previously been proposed [35,36], as well as the start of rev (and hence a section with overlapping open reading frames and resultant conserved wobble positions) and a GAR exon splicing enhancer [37,38]. Given the proximity of the regions 6206–6226 and 6324–6611, we folded the entire region 6206–6611 to look for correspondences between elements within the larger region (S15 Fig). This fold gives a correspondence between the regions 6218–6229 and 6597–6608, with a helical structure in a long-range interaction. Stem-loop structures have been postulated 3′ to the env signal peptide [2,39], and we considered whether there may be splicing-dependent folding with the potential to affect translation and/or function. To this end we have produced local folds of alignments corresponding to the three regions of interest with a region directly 5′ to them (S16 Fig), and of alignments corresponding to spliced RNA of the 5′ leader region with the donor site D1 spliced to each of the splice acceptor sites A4c, A4a, A4b, A5, A5a and A5b, and the remainder of the region HXB2 5951–6611 lying 3′ to the respective splice acceptor (S17–S22 Figs; see S22 Table for additional information on the sequences used to produce these figures). These folds indicate that the region HXB2 6206–6226 may have the ability to form stable plausibly functional alternative structures by base pairing intramolecularly in the unspliced RNA or in the chimeric individual spliced species that arise from use of the important A4/5 splice site cluster (Fig 7). Our folds also suggest that the region HXB2 6324–6611 contains a number of conserved structural motifs. Our method highlights a number of further regions of interest, some corresponding to known features of the HIV-1 genome but others whose conservation we cannot fully explain. We detail these below in the sequence in which they appear in the genome, for ease of reference. The output of nMPD analysis of A-subtype sequences is shown in S28 Fig. S21 Table summarizes regions of interest found and S29 Fig. The validation analysis reproduces all of the regions found in the main analysis with the exception of one region of pol, which can be seen in the plot of cumulative analysis in S29 Fig to be conserved, but is not deemed significant by the algorithm. Conversely, the validation analysis finds two regions in gag that it deems to be significantly conserved, which can be seen to be conserved in the main B-subtype analysis (S10 Fig), but are not deemed significant in the main analysis. We also briefly analysed the coding potential of the possible alternative ORFs we have described in pol and env. For the possible alternative ORF in pol, 95 of the 107 sequences that we studied (89%) contained appropriate initiation codons; for the possible ORF in env, 106 of the 107 sequences that we studied (99%) contained appropriate initiation codons; the spread of stop codons was similar to that in the B-subtype sequences. Again, there is consistent conservation of the appropriate Kozak consensus nucleotides for the possible open reading frame, with high purine conservation at the −3 positions and guanine conservation at the + 4 positions relative to position 1 of the AUG codons (S30 Fig).
We have applied to a large collection of collated, previously experimentally acquired HIV-1 genomic data a novel version of an informatic pipeline for detecting conserved genetic regions. This combines a method for detecting nucleotide conservation independently of amino acid constraints with a method for detecting runs of conservation without requiring prior knowledge of the scale of any feature of interest. Our pipeline reassuringly detects features where it is already known that nucleotide conservation is required—for example, the Rev-response element is readily identified, and there is a striking correlation between regions our pipeline identifies and known clusters of splice sites—but it goes on to detect novel features whose role has not been described previously. By combining this pipeline with structural prediction tools capable of accounting for sequence variation in the prediction made, we propose a series of conserved structures, some with similarity to structures already published, but some that are novel and plausible, including a highly stable structure surrounding the central polypurine tract, which parallels a structure surrounding the 3′ polypurine tract. The structure we show (Fig 6) would permit exposure of the central polypurine tract for priming of second strand synthesis [62,63] when it forms a DNA/RNA hybrid: previous work proposes that an oligo-A structure such as this in the context of an inverted repeat would promote reverse transcriptase pausing or even dissociation from the template [64,65]. Pausing would allow time for second strand synthesis to start using the central polypurine tract RNA as a template before it is hydrolysed; dissociation may leave the central polypurine tract RNA in situ for use as a template. Formation of G tetrads could provide a switch mechanism to an alternative conformation that would disrupt the structure we propose when the RNA is in the monomeric or dimeric form, and might have roles in obscuring double-stranded RNA runs to evade innate immunity or to foster intergenome dimer formation for encapsidation. In one case (the tat/vpu/env region), the structures we propose, taken in conjunction with the identified regions of nucleotide conservation, immediately lend themselves to plausible functions in alternative splicing. Variable bridging of the HXB2 6206–6226 sequence between the 5′ leader region or the unspliced region upstream of HXB2 6206 might either promote translation of the partially spliced env transcripts from the immediately downstream env translational initiation codon, or stabilise part of the structure located 5′ of this region that contains the splice donor site D4, which is used in conjunction with downstream splice acceptor sites to generate completely spliced mRNA transcripts. Additionally, the positioning of the conserved regions of interest we have identified within the spliced folds, together with the folds themselves, raises the possibility that splicing from the D1 donor site to the acceptor sites in the A4c to A5b region may result in a conformational stabilisation of the 5′ cap, possibly contributing to the extensive utilisation of the A4/5 splice acceptor cluster for over 90% of env, vpu, rev and nef transcripts [66]. The leader region can adopt two alternating structures [67], and stabilisation of the leader region in the long distance interaction (LDI) conformation, rather than the branched multiple hairpins (BMH) conformation [68], may imply a function for this elongated helix loop in the spliced RNA. The LDI conformation is proposed to reduce inter-molecular genome dimerisation and this would be a valuable property to preserve in spliced mRNA to prevent dimerisation using the high-affinity SL1 dimer linkage site, restricting this property to unspliced viral RNA. Although a previous in vivo analysis of RNA structure in this region did not identify the LDI conformation [69], detection of this in spliced RNA may have been compromised by use of a primer downstream of the gag start site. The key to generating these structures—and the novel step we have undertaken—is priming the folding algorithm with just the conserved region, identified by our scale-free analysis of nMPDs, that may form a conserved structure, to increase the algorithm’s ability to select and predict the locally conserved secondary structure. To prime in this way cannot be achieved without knowledge of the extent of the conserved region, which is a novel feature of the detection pipeline we have used. Whilst we have strong mathematical evidence for the existence of conserved regions we rely on prior publications to support the structural predictions we have made, where such evidence exists. The largest and most striking region of conservation we identify is the Rev-response element and this is powerful validation of the predictive accuracy of our modelling. Other regions we show are based on what are highly accurate structural prediction programs but some need functional and structural validation and must be accepted as putative at present. We have mentioned two regions where one explanation of genetic conservation would be hithertofore undescribed open reading frames. We caution that we are not aware of any empirical evidence of translation of these putative frames. The lack of strong database search hits against peptides of known function indicates that putative products of these frames are not directly analagous to already-characterised peptides. Some regions of the genome where ovelapping reading frames are known to exist are not identified as having a high degree of conservation by our method, indicating that the degree of variability of some regions of HIV-1 without wobble positions is comparable with the background degree of variability. Our method relies upon what is in effect a signal processing algorithm that most easily picks up relatively long runs of conserved genetic elements in a background where there is less conservation. The algorithm is applied gene-by-gene. This explains why, for example, the algorithm does well at detecting conserved elements within the env gene, where the overall high degree of variability makes it easier to identify conserved elements (in signal processing terms, the signal-to-noise ratio is high). It also explains why, for example, the algorithm fails to detect the well-characterised frameshift structure towards the end of gag, which has only a relatively short run of highly conserved nucleotides. It is important to understand that this algorithm is optimised for identifying longer contiguous conserved genetic elements such as those where RNA structure needs to be conserved and it may miss short or non-contiguous motifs, especially in genes with relatively lower nucleic acid variability overall. Whilst there is no minimum run length the algorithm can detect (it can pick up any sufficiently strong signal of length at least 2), the degree of conservation of a region versus background variability needs to be much higher for detection of shorter conserved regions. This is illustrated in Fig 5, left panel of our paper describing the algorithm [13], and to aid understanding of the degree of conservation required to detect much shorter regions using this algorithm, we have extended our benchmarking simulation to illustrate results for signal lengths down to 2 (S31 Fig). Our validation step using A-subtype sequences shows consistency of the method’s performance between related datasets, whilst illustrating a feature of our algorithm that minor differences between datasets may lead to minor differences between the start and end positions of regions deemed significantly conserved. The nMPD algorithm is not applied to non-coding regions at all and will not pick up conserved elements in non-coding regions. Thus overall our method cannot detect known microRNAs or long non-coding RNAs [70]. This novel approach to detection of nucleotide conservation in HIV-1 has revealed unexpected insights into the viral genome and can be applied similarly to other RNA viral sequences. The identification of unsuspected regions of conservation is of importance as it directs experimental approaches for investigation into regions of conserved structure or function. In conjunction with further experimental validation, it can suggest loci that are likely to represent stable drug targets.
HIV-1 is a very rapidly mutating organism, however some parts of its genetic material change more than others. We looked for coding regions of HIV-1 that change relatively little, by turning the problem of finding such regions into a problem in signal processing, and solving this using a novel analytical approach that we recently described. We investigated why the regions we identified change less, including using the genetic code in the regions we found to prime an algorithm to predict their structures. In some cases there are already known functions for the features we found, in others they provide new insights into the properties of known regions, and in some cases we identify new regions that vary less for as yet unknown functional reasons.
sequencing techniques
medicine and health sciences
pathology and laboratory medicine
pathogens
microbiology
retroviruses
viruses
immunodeficiency viruses
rna viruses
genome analysis
sequence motif analysis
molecular biology techniques
research and analysis methods
sequence analysis
rna structure
genomics
lentivirus
sequence alignment
bioinformatics
medical microbiology
hiv
microbial pathogens
hiv-1
molecular biology
nucleotide sequencing
biochemistry
rna
rna folding
dna sequence analysis
nucleic acids
database and informatics methods
viral pathogens
genetics
biology and life sciences
computational biology
gene prediction
organisms
macromolecular structure analysis
A scale-free analysis of the HIV-1 genome demonstrates multiple conserved regions of structural and functional importance