We chose HIV-1 group M as our model system of study because the group M subtype is the most predominant HIV-1 strains which are responsible for AIDS pandemic. HIV-1 group M reference sequences were downloaded from http://www.hiv.lanl.gov. The list of HIV-1 group M subtype genomes used in this study is provided in Supplementary Table 1. All of these group M subtype genomes were of non-recombinant types. The genomic coordinates of each gene from each of the HIV-1 genomes were downloaded from GenBank. From these genomic coordinates, we have identified regions of gene to gene overlap in the respective genomes (Sabath et al., 2008). We considered those regions as true overlaps where a two adjacent genes involved in overlapping are in different reading frames. The details of overlapping gene pairs of each HIV-1 genome used in this study, along with their genomic coordinates and reading frames are enlisted in Supplementary Table 1.

The residue wise disorder score for each protein was predicted by IUPred algorithm (Dosztanyi et al., 2005a) using the option to predict long disordered regions. We choose IUPred algorithm for disorder prediction as IUPred didn’t train on any specific dataset and hence it gives an unbiased estimate of disorder score (Dosztanyi et al., 2007). Moreover, this algorithm has been validated to reflect truly the structural status of proteins in viral genomes (Hagai et al., 2014). Based on the pairwise interaction energy score, IUPred assigns disorder scores for each amino acid (Dosztanyi et al., 2005b). Residues having a disorder score above threshold value 0.5 were considered as disordered. From these predicted disorder scores, we calculated several measures of protein disorder content. We have measured the percentage of disordered residues by calculating the total disordered residues to the total amino acid residues of the proteins. Next, the percentages of disordered residues were averaged for all proteins in the proteome.

The short linear motif patterns were downloaded from ELMdb¹ (Dinkel et al., 2012). The prediction of SLiMs on HIV-1 proteome was done using in-house perl scripts. In order to predict true SLiM in OVR as well as in non-overlapping region (NOVR) we took an approach similar to Hagai et al. (2014) where, we have shuffled the entire protein sequence of each gene 1000 times and predicted SLiMs in the shuffled set of sequences. The SLiMs that were retained in both the shuffled set of sequences as well as the original type were referred to be false positives, while the SLiMs that occurred in original sequences but disappeared upon shuffling were referred to as the original or true SLiMs. We scored each SLiM according to their occurrences in the randomly shuffled sequence set, for instance if a SLiM occur 10 times out of 1000, then the probability of occurrence of that SLiM in shuffled set is 0.01. Next, we classified those SLiMs as true SLiMs which has a probability of occurrence in the shuffled set is less than 0.1. The sequences were shuffled using the Ushuffle program (Jiang et al., 2008) with seed 10^-4. Experimentally validated SLiMs that are involved in interaction with host proteins were retrieved from the published dataset of Davey et al. (2011) and Halehalli and Nagarajaram (2015).

We calculated SLiMs which are present within intrinsically disordered residues, since they may act as potential molecular recognition features (Mohan et al., 2006). We identified true SLiMs whose positions coincided with intrinsically disordered residues and termed them as functionally potential SLiMs.

Conservation levels of the respective true SLiM positions were estimated where a conservation score was assigned to each of the SLiMs identified. The conservation score was calculated as the ratio of the number of HIV-1 genome in which a given SLiM was retained in the same position in all orthologous genes of the respective genomes and the total number of HIV-1 genomes studied.

Prediction of protein phosphorylation (PP) sites were done using NetPhos 2.0 (Blom et al., 1999) which gives accurate predictions using neural network for serine, threonine and tyrosine residues in a given protein. It assigns a score to each serine, threonine or tyrosine residues indicating the probability of these residues of getting phosphorylated. Residues with scores greater than 0.50, indicate potential PP sites. List of experimentally validated phosphorylation sites were obtained from (Schwartz and Church, 2010).

Mann–Whitney U test was used for pairwise comparison of different factors used in this study. Spearman’s bivariate correlations were used to find out the association between factors. All of these statistical analyses were done using SPSS 13.0 package.

Previously, it was reported that the HIV-1 genome contains overlapping genesthat contribute to synonymous site conservation in the virus (Mayrose et al., 2013). Here, we identified 235 OVRs of varying lengths across 35 HIV-1 group M subtype genomes (Supplementary Table 1). The proportion of OVRs of a given genome was then calculated by summing up all the lengths of OVRs present in that genome, and then dividing it by the total length of the coding sequence of that genome. We observed that proportion of overlap varies from 4.8 to 8.0% of the total coding regions (Figure 1; Supplementary Table 1). The density plot (Figure 1) demonstrates that majority of the genomes have proportion of OVR of 8% of the total coding regions.

It was reported that viruses of eukaryotes, often encode for proteins with unstructured regions those are less constrained and more likely to accumulate mutations (Pushker et al., 2013). This may in turn expedite the advent of novel functions during the evolution of the virus (Gitlin et al., 2014). Previously, it was shown that OVR in unspliced RNA viruses encode an elevated degree of protein intrinsic disorder (PID; Rancurel et al., 2009). So, we intended to investigate whether the OVRs of HIV-1 genome have the potency to encode PID. In calculating the proportion of disordered residues within the OVRs and NOVRs of all the genes of 35 HIV-1 group M subtype genomes, it was noticed that disordered residues are significantly (P = 0.003) overrepresented in OVR than NOVR (Figure 2A). To study the effect of gene overlapping on PID more explicitly, we decided to split the data into eight HIV-1 genes (tat, rev, gag, pol, vif, vpr, env, and vpu). It was observed that except for vpu (which has a very ordered and stable secondary structure and a very little number of disordered residues in OVR), all other proteins encode considerable amount of PID in OVR as compared to that in NOVR (Figure 2B). In addition, we calculated the proportion of intrinsically disordered residues for each gene of 35 HIV-1 subtype genomes. It was observed that the proportion of OVRs of each gene shared a strong, significant positive correlation with the proportion of predicted intrinsically disordered residues in the protein encoded by the corresponding genes (Figure 2C).

Moreover, OVR stretches are short and terminal in nature, and short sequences are predicted to be highly disordered. We next decided to investigate whether the increase of structural disorder in OVR is independent of the length of the protein. We did a linear regression analysis taking total disordered residues in a given protein and its length as independent covariates and OVR length as dependent variable. The result of the multiple linear regression analysis is represented in Table 1. It could be interpreted from the given multivariate linear regression that increase of intrinsic disorder in OVR is independent of protein length.

Next, we planned to validate whether regions predicted by IUPRED is actually disordered or not. For this purpose we retrieved the structure of individual HIV-1 proteins from Protein Databank (PDB) and tried to ascertain whether the OVR in individual proteins are actually disordered. It was previously shown by Le Gall et al. (2007) that the missing residues (unobserved residues) in X ray crystallographic structures could be structurally disordered, due to their highly dynamic atom positions which vary significantly over time. Taking this into consideration, we have mapped the OVR in PDB structure of HIV-1 proteins, the result of which are delineated in Supplementary Table 2. It was observed that OVRs are overrepresented in the “unobserved” or missing residue category of PDB structure. Thus, it is reasonable to assume that OVR residues could be highly dynamic and its structure could not be deciphered using standard X-Ray crystallographic protocols. Moreover, the structure of vpr and vpu are delineated in Figures 2D,E with their OVR highlighted in red. It could be observed that the highlighted region, i.e., OVR are structurally more flexible.

In the previous section, we observed that OVR encodes significant amount of PID. It has been previously reported that PID plays pivotal role in protein–protein interactions (Uversky, 2013). Therefore, we were curious to investigate whether OVRs play any role in virus–host interaction or not. In the recent years, there are considerable volume of researches that discovered a profound role of SLiMs in viral proteome that readily facilitate virus host protein–protein interaction (Hagai et al., 2014). Thus, preponderance of these types of motifs in a given viral sequence often acts as an indicator of its involvement in protein–protein interaction. Henceforth, we tried to ascertain whether OVRs in HIV-1 genome contains SLiMs to initiate virus–host interaction. In this context, we intended to find out whether OVRs, apart from encoding structural disorder are also involved in encoding SLiMs that could be essential for host–viral interactions. We scanned all the genes involved in overlapping across 35 HIV-1 group M genomes and figured out the position of different SLiMs along the entire protein length. Similarly, we shuffled individual protein sequences and figured out the probability of occurrence of SLiMs in the shuffled set. The SLiMs that occurred at a low frequency or does not occur at all in shuffled set were considered as true SLiMs. It was observed that OVRs are enriched in true SLiMs as compared to NOVR which has a depleted density of true SLiMs. (Figure 3A). Moreover, it was observed that the proportion of OVRs of a gene shares a significantly positive correlation with the SLiM density of the concerned genes (Figure 3B). We next asked whether SLiMs residing in OVRs were functionally more important than the ones residing in NOVRs. Previously, Hagai et al. (2014) determined the functional potency of SLiMs by calculating the number of SLiMs residing in disordered region of proteins. They have proposed that SLiMs, residing within locally flexible disordered regions are more suitable to participate in protein–protein interactions and thus these SLiMs are termed as functionally potential SLiMs. Similar to their analysis, we also computed densities of functionally potential SLiMs (SLiMs localizing within structurally disordered residues of proteins) in OVR as well as in NOVRs of HIV-1 genomes. It was observed that densities of functionally potential SLiMs were far higher in the OVR as compared to the NOVR (Figure 3C).

We further compiled a dataset of experimentally validated interaction SLiMs that are reported to interact with host proteome. We retrieved the list of genes and their corresponding interaction SLiMs Supplementary Table 3a), and found out whether these SLiMs are located within OVR or NOVR of that gene. Out of 71 SLiM-gene pair only 21 unique SLiM-gene pairs were observed to be matching with our true SLiM dataset. It was observed that interacting SLiMs were enriched in OVR for tat and rev proteins, gag protein also contained a few experimentally validated SLiMs in their OVR. In contrast, vif, vpu and vpr, and env have all their experimentally validated SLiMs in NOVR (Supplementary Table 3b).

Next we intended to find out whether the SLiMs discovered in the OVR of HIV-1 proteins are well conserved in their position across all subtypes of HIV-1 genome. For this, we have estimated the conservation level of each SLiMs located within OVRs as well as in the NOVRs. It was observed that there was no significant difference of conservation level between SLiMs located in the OVR to that in the NOVR Figure 3D. This may be due to higher intrinsic disorder content of OVR which relaxes the level of purifying selection on these SLiMs.

Another potential indicator of virus–host interaction is PTMs; Hundt et al., 2013). Precisely, among the different PTMs that take place in vivo, PP is one of the most important and well studied types of PTM. In previous works, PP has been shown to play a pivotal role in HIV-1 interactome (Francis et al., 2011). Therefore, we were interested in exploring the role of gene overlapping in effecting PP across the protein length. In order to ascertain this we have mapped the OVR of each gene and calculated the PP density of OVRs as well NOVRs. It was observed that the PP densities of OVRs were far greater than that of NOVRs (Figure 4A). It was also observed that the proportion of gene overlapping correlated with a fraction of PP density of a given gene (Figure 4B). We have confirmed our results using experimentally validated PP sites, the results of which are delineated in Supplementary Table 4.

We have found that OVRs attributed to the overrepresentation of SLiMs and PPs in HIV-1 genome. Now, it was also reported that SLiMs and PP sites encompass affluence of disordered residues (Iakoucheva et al., 2004; Meszaros et al., 2012). Thus, we have asked whether the enrichment of SLiMs and PP sites are by-products of an abundance of intrinsic disorder in the respective OVRs or OVRs by virtue harbors more SLiMs and PP sites as compared to NOVRs. In other words, we were interested to find out whether the overrepresentation of SLiMs and PPs in OVRs are independent of PID. For this, divided each HIV-1 subtype proteome into three categories- (i) OVRs with 80% or more of its residues in disordered state (DO) (ii) OVRs with 80% of its residue in the structured state (SO) (iii) NOVR with 80% or more of its residue in structured state PID (SN). Our hypothesis, in this context was that DO would entail elevated degree of SLiMs and PPs as compared to SO and SN, since it is highly disordered as well as overlapping in nature. In addition, SO would entail an elevated density of PP and SLiM as compared to SN if the occurrences of PPs and SLiMs are mutually exclusive of PID. This is because the phenomenon of gene overlapping will itself be responsible for enrichment SLiMs and PPs in SO, whereas SN (where there is no genetic overlap) would definitely confer depleted density of SLiMs and PPs. On the other hand, the absence of significant difference of PP and SLiM density between SO and SN would reflect that PID is the reason behind the elevated densities of PPs and SLiMs in the OVR and OVR doesn’t have any significant impact over PP and SLiM density. Finally, we estimated the SLiM and PP density for these three categories separately, the result of which are delineated in (Figure 5). It could be observed that SLiM content is higher in SO as compared to SN which reflects that acquisition of elevated frequency of SLiMs in the OVRs is mutually exclusive of PID. On the contrary occurrence of PP was not mutually exclusive of PID since there was no significant difference of PP density between SO and SN groups. Hence, it could be interpreted that although acquisition of SLiMs in OVR is exclusive of PID, but enrichment of PP in OVR was mainly due to increased disorder content of OVRs.