1.
Introduction
Computational methods can play an important role in bridging the sequence–function–structure gap [1]–[3]. With a series of structure prediction programs helical structures within globular proteins [4] or structures such as helical transmembrane domains (TMDs) [5],[6] can be predicted. Far fewer programs are available predicting β-sheet as the membrane spanning motif. Beyond this level of predicting the secondary structure, obtaining the tertiary or even quaternary structure is within reach [7]. In terms of obtaining functional information, sequence alignment is a sophisticated tool to compare unknown sequences with those, for which for example more information upon its function and structure is available [8],[9]. The alignment can be done solely on the sequence as such, or focused on particular sequence motifs, such as charged amino acids in DNA binding proteins, to obtain spatial information about the protein [10].
The combined search, for secondary structure as well as functional information, is applied to the viral non-structural protein 4 (NS4B) of hepatitis C virus (HCV). NS4B is a 261 amino acid membrane protein. It is expressed in the infected cell as part of the approximately 3000 amino acid long polyprotein [11]. Consequent cleavage of the polyprotein by viral and cellular proteases release the protein together with 9 others. The role of the polytopic membrane protein is to support the generation of viral replication organelles [12]. NS4B is proposed to harbor two amphipathic helices on either side of its core region. The core region harbors at least four TMDs based on secondary structure prediction programs (SSPPs) applied to this region [13]–[15]. TMDs 2 and 3 are suggested to be held together by a nucleoside-triphosphatase motif [16],[17]. Eventually a fifth TMD is also proposed from experimental results [18]. The oligomeric state is experimentally proposed to be that of at least a dimer or trimer with higher oligomeric states [19].
A sequence of steps is proposed to make a prediction of the putative structural and functional features of NS4B by applying bioinformatics tools. Secondary structure prediction programs (SSPPs) are used to identify the putative structural motif of the TMDs and their length. The proposed sequences of NS4B is aligned with proteins available in the protein structure and sequence data base (Protein Structure Data Bank, UniProt) to identify reasonable sequence overlap and to propose functional features. In addition, with the outcome of the alignment the search is narrowed into the search for functionally related proteins. This also includes membrane protein Nsp3 of corona viruses. Together with the functional and structural analysis the functional role of individual TMDs is suggested. The data are compared with those obtained when using the reported sequence of TMDs (hither forth denoted as r-TMDs) in the literature [17]
2.
Materials and methods
2.1. Secondary structure prediction
The sequence of NS4B is based on that of HCV strain JFH-1(2a), GenBank accession number AB047639 [20] and used to predict consensus TMDs, hither forth called p-TMDs, using SSPPs such as CCTOP[21], DAS (cutoff@1.7) [22] and DAS-TMfilter [23], HMMTop [24], MemBrain [25], Memsat [26] and MEMSAT3 [27], OCTOPUS [28]–[30] and SPOCTOPUS [28],[29], Philius [31], Phobius [28] and PolyPhobius [32], Pro and Prodiv [33], Scampi [34] and ScampiMsa [34], SPLIT 4.0 [35], TMHMM 2.0 [36], TMMOD [37], TMpred [38], TOPOCONS [39], PSIPRED 4.0 [40].
The reported sequences of the TMDs reported in the literature (r-TMDs) were as follows [17]: r-TMD1 70–90: LPGNPAVASMMAFSAALTSPL; r-TMD2 112–131: PAGATGFVVSGLVGAAVGSI; r-TMD3 137–154: LVDILAGYGAGISGALVA; r-TMD4 172–190: LPGILSPGALVVGVICAAI.
Targeted TMDs used for alignment were taken from membrane proteins involved in Ca-homeostasis and grouped as (i) Ca-ion conducting channels: the calcium release-activated calcium channel (Orai) from Drosophila melanogaster (PDB ID: 4HKR), L-type, high voltage-activated Cav channel (Cav1.1) from rabbit (PDB ID: 3JBR), human voltage-dependent L-type calcium channel (Cav1.2) (Uniprot Q13936), human ryanodine receptor 1 (RyR1) (Uniprot P21817), human sarcoplasmic/endoplasmic reticulum Ca2+-permeable cation channel inositol-1,4,5-triphosphate receptor type 1 (IP3R1) (Uniprot Q14643); (ii) Ca-transporters: mitochondrial calcium uniporter (MCU) from C. elegans (PDB ID: 5ID3), rabbit skeletal-muscle sarcoplasmic/endoplasmic reticulum Ca2+-ATPase (SERCA) (PDB ID: 3TLM); (iii) Ca-mediators: SERCA inhibiting rabbit phospholamban (PLN) (PDB ID: 2KYV), human stromal interaction molecule 1 (STIM1) mediating store-operated Ca2+ entry (Uniprot Q13586); and (iv) Nsp3 of coronaviruses being involved in the formation of replication vesicles: mouse hepatitis virus (MHV) (Uniprot: Q66WN6, a.a. 2062-2841), SARS-CoV (Uniprot: P0C6X7, a.a. 1997-2745), SARS-CoV-2 (NCBI Reference Sequence: YP_009725299.1, a.a. 1202-1945), MERS-CoV (Uniprot: K0BWD0, a.a. 1953-2750).
2.2. Sequence alignment and structure prediction
Protein-BLAST (https://blast.ncbi.nlm.nih.gov) was used to find similar sequence segments from existing structures in its default mode. The searching was based on the default database excluding data bases like Models (XM/XP), Non-Redundant RefSeq Proteins (WP) and Uncultured/Environmental Sample Sequences. The maximum target sequences were set to 20,000. Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/) was used in its default mode to search for proteins with similar function by sequence alignment and to create the phylogram. Sequence Identity and similarity of each aligned sequence pair was calculated using Software Sequence Manipulation Suite (SMS) (http://www.bioinformatics.org/sms2/index.html). All programs were used in their default mode.
3.
Results
An algorithm is proposed to predict structure – function relations of the transmembrane domains (TMDs) of the viral membrane protein NS4B (Figure 1). Taking the primary structure of the protein as input from outside sources, e.g. experiments, proposals of the stretch of putative secondary structural motifs is made. With the stretch of amino acids including the secondary structural motif a two-level sequence alignment is applied. On a first level, proposal of the function of the investigated protein is made by comparison with proteins of known function which, at best, also includes structural information. On a second level, with the knowledge from the first level, the search space is narrowed to focus only on those proteins with specific functions and in this study on their TMDs. Other information from outside experimental sources can also enter at this level. The algorithm is not limited to the protein investigated in this study.
3.1. Secondary structure prediction
A series of 20 SSPPs is used to search for a consensus region of the putative TMDs within the sequence of NS4B (Figure 2). A consensus is reached for a particular amino acid to be part of the TMD if more than 19 of the SSPPs propose this amino acid to be in a helical conformation (see also Figure 3 for details). When using the full-length sequence of NS4B consisting of 261 amino acids, the best consensus is reached for a stretch around I169 to H194, here denoted as p-TMD4 (e.g. red box in Figure 2). A stretch from A110 to G160 is also predicted by all of the programs to containing a helical motif. By visible inspection of the chart (Figure 1, I), it is intriguing to define a short loop for G132 to D139. Thus, two TMDs are predicted for this region, denoted as p-TMD2 and 3 (e.g. red boxes in Figure 2). This observation is supported by the fact that the full stretch with its 50 amino acids is too long to be a TMD, which is usually in the range of 17–25 amino acids [41],[42]. Another helical motif is proposed between residues W-50 to Q-108 by 15 of the SSPPs (Figure 2,I). From those, two of them even propose two helical domains. The stretch overlaps with the stretch proposed to harbor AH2 (S42–G66 [17]). By visible inspection, a consensus sequence M80 to Q108 is chosen as TMD1 which includes the prediction of most of the programs as well as the putative second helix (e.g. red box in Figure 2).
Using the individual sequence which potentially harbors p-TMDs 2 and 3, A106 to M165, supports the existence of two TMDs as mentioned above (Figure 2, II). Using a separate stretch V46 to A110 enhances the impression of the existence of two possible TMDs.
Focusing on the individual stretches harboring potentially the TMDs reveals that the prediction for p-TMDs 2 (blue), 3 (green), and 4 (purple) are suggested for most of the programs within a narrow region of the amino acid sequence (Figure 2, III). The stretch for p-TMD1 is also predicted when including the amino acid sequence from I61 to A110 (Figure 2, III, red). Using a separate stretch A41 to M80, nine programs propose the aforementioned second helical stretch (Figure 2, III, orange).
A similar study using nAChR, GLIC and LMP1 of Epstein-Barr virus shows that independent of the length of the amino acid sequence used, the prediction of the TMDs is similar (Suppl. Figue 1 and 2). In most cases the predictions coincide with those proposals found based on experimental evidence in the literature (black boxes in Supplemental Figures 1 A,B,C).
Counting all the hits for the amino acids being in a helical conformation from the 20 programs defines p-TMD4 of NS4B from L171 to L191 (Figure 3). This definition and those for the other TMDs are shown as red boxes in Figure 1. Also, the region containing p-TMDs 2 and 3, P112 to G160, show distinct patterns for two maxima as well as the N terminal region (W50 to I109). Allowing for a threshold that at least 11 programs suggest a helical stretch the following TMD stretches are proposed:
p-TMD1 from F82 to G102 (FSAALTSPLS TSTTILLNIM G),
p-TMD2 from G117 to L133 (GFVVSGLVGA AVGSIGL), and
p-TMD3 from I140 to F155 (ILAGYGAGIS GALVAF), and
p-TMD4 from L171 to L191 (LLPGILSPGA LVVGVICAAI L).
Residues G134 to D139 (GKVLVD) are seen as a short loop between p-TMD2 and 3 based on their potency to form salt bridges. When using the sequence of the three loops between the p-TMDs (G103–T116, G134–D139, K156–N170) neither sequence alignment (BLAST) nor motif search servers (InterPro [43], PROSITE [44]) identify any similarity with membrane proteins.
The percentage of overlaps of the p-TMDs with the r-TMDs reported in the literature [17] (see black boxes in Figure 2) decreases in the order TMD4 (90.5 %) > TMD3 (79 %) > TMD2 (68.2 %) > TMD1 (27.4 %).
As a result, while prediction of TMD2, 3, and 4 are fairly similar, there is some ambiguity for predicting TMD1.
3.2. Sequence alignment and structural Bioinformatics with experimentally derived sequence data
The software BLAST is used to identify proteins which show sequence identity with especially the TMDs of NS4B and for which structural information is available. With the latter information a structural functional proposal of the p-TMDs of NS4B is given. Using full length sequence, NS4B [20] aligns with 80 % identity with itself when embedded in a polyprotein (P27958) (Table 1,i).
Limiting the subject sequence to the p-TMDs containing segments, F82–L191 of NS4B, results in no alignment (Table 1,ii). For the stretch including the r-TMDs, L70 to I190, an enzyme, succinate-CoA ligase, is identified (Suppl. Table 1).
Aligning individual ranges which include the individual sequences of the p-TMDs, identifies sequence identities with membrane proteins in the range of 67–100 % identity (Table 1,iii). The top scored protein is a chromate efflux transporter with an identity of 100 % for p-TMD2, followed by an ABC transporter permease (89 %) for TMD3. All other membrane proteins are mainly transporters for the other p-TMDs. p-TMDs 1 and 4 also show sequence identity with a transferase (71 %) and ion channels (Vpu of HIV-1 (71 %) and an unknown ion channel (67 %)), respectively. Using the individual r-TMDs, transport proteins are proposed for r-TMD2 and 4 as well (Suppl. Figure 1). Selecting r-TMDs2 and 4, for the latter also a Ca/Na antiporter is identified.
As a result, when focusing on the TMD sequences, an increased overlap with the TMDs of other membrane proteins such as transporters, transferase and ion channels is achieved. One of the overlapping proteins, is a Na+/Ca2+ exchanger, involved in Ca-homeostasis. This supports evidence from experiments, showing that NS4B could be involved in Ca-homeostasis of the endoplasmic reticulum [45].
3.3. Sequence alignment with function-related structure-based sequences
Three classes of membrane proteins, (i) to (iii) (see Materials and Methods) knowingly involved in Ca-homeostasis are chosen, to probe functional correlation of their reported functional TMDs (Figure 4) with the p-TMDs of NS4B. Values such as identity (I) and similarity (S) are obtained for the p-TMDs and the selected proteins. The values for each of the p-TMDs for all the 11 proteins are plotted as single values or as an averaged if the protein consists of a few subunits to identify the differences between the TMDs. The values for I are in the range of 5.1 (Orai) to 20.2 ± 5.9 (IP3R), those of S in the range 14.0 (Orai) to 52.2 (PLN) (Table 2).
Counting the matches of the p-TMDs of NS4B which have overlap with either pore lining TMDs or are identified to be involved in ion translocation at the protein site, focusing on classes (i) and (ii) proteins identifies the sequence as follows. There are 3 matches for p-TMD1, 4 matches for p-TMD2, 3 matches for p-TMD3 and 8 matches for p-TMD4 when looking at I and S values (Table 2). p-TMD2 and 4 overlap with TMDs of PLN as well as STIM1, respectively. A similar distribution is also obtained for the r-TMDs (1 match (r-TMD1), 6 matches (r-TMD2), 1 match (r-TMD3), 8 matches (r-TMD4) (Suppl. Table 2). For classes (iii) r-TMD2 and 3 have overlap with domain II of PLN, and r-TMD4 with STIM1.
Vpu of human immunodeficiency virus type 1 (HIV-1) is identified as a member of the class of channel forming proteins, also called viroporins, with low ion selectivity [46],[47]. Analysis of the sequence alignment on the basis of the individual p-TMDs with specifically the TMD of Vpu (Figure 5) reveals that p-TMD4 has the largest I value of 25.9 % followed by p-TMD1 with 22.2 % (Suppl. Table 3). On the basis of S values, p-TMD2 with 56.5 % stands out compared to p-TMD4 of 51.9 %.
TMDs 2 and 4 are proposed to have an active role in the function of NS4B. Small changes in the sequence, as for p-TMDs versus r-TMDs, do not affect the results very much. It is like a second experiment under similar conditions.
3.4. Sequence alignment with Nsp3 of coronaviruses
The sequences of the two TMDs in each of the four Nsp3 proteins are all predicted within narrow stretches of the amino acid sequence by the secondary structure prediction programs (Suppl. Figure 4). Since the prediction of TMD2 includes parts of the proposed ECTO domain [48], two type of stretches are chosen for pairwise alignment, one which refers to the stretch proposed by the secondary structure prediction programs which then overlaps with the stretch found to be part of the ECTO domain (TMD2*), and one stretch in which this part is not included (TMD2).
The sequences for both TMDs hypothesizes a sister group for SARS and SARS2 with differences in sequence similarity of 12–19 % (Suppl. Figure 5). MERS is shown as an outgroups sequence difference of 30–39 % to the sister group especially for TMD2* and TMD1, respectively (Suppl. Figure 5A,B). From this clade, MHV is hypothesized to have 36–40 % differences. For TMD2-(-E) the SARS sequences and MHV for a clade (33 %) from which MERS shows a difference of 40 % in the sequence (Suppl. Figure 5C).
Individual alignment of the four TMDs of NS4B with the two TMDs of Nsp3 of coronaviruses (Clustal Omega) (for sequence see Suppl. Tab. 4) reveals that p-TMDs 2 (7 times) and 4 (5 times) of NS4B show the largest number of highest I / S values amongst its four TMDs (Suppl. Table 6, and Suppl. Table 5). Using the stretch which includes the ECTO domain, p-TMD1 (4 times) shows the largest increase in number of hits (+3) compared to the value without the ECTO domain. It also shows highest numbers together with p-TMD 2 (4 times) and TMD4 (5 times). Choosing the r-TMDs for alignment (no ECTO domain) TMD1 (7 times) and r-TMD2 (6 times) are most often the highest values for I and S (Suppl. Tables 7 and 8). Including the ECTO domain reveals that r-TMDs 1 and 3 (6 and 5 times, respectively) show the highest number of hits with the latter showing the largest increase (+3) amongst the TMDs of NS4B.
p-TMDs 1 and potentially 3 show similarity to the ECTO domain of Nsp3. p-TMDs 4 and 2 are likely to be indifferent to the changes in the alignment protocol.
4.
Discussion
With the advance of sequencing the genome rapidly an enormous amount of protein sequences is at hand. Increasingly slower is the pace of converting this information into functional, then structural and even dynamical information of the protein leading to the sequence–structure–function gap [2]. This gap is even more pronounced for membrane proteins, despite their importance in current therapy [49]. Computational tools can be seen as a complementary tool to the experimental techniques to enhance the speed of delivery of essential information [1]–[3]. Besides generating software able to predict particular features it is also essential to develop algorithms which can be used to achieve meaningful information. In this respect, a flow-chart is presented using bioinformatic tools to achieve potential structure–function correlation. This chart predicts the structural features of the target protein on a first level, here e.g. transmembrane helical motifs, and correlates this motif with information available in sequence and structure data banks on the basis of protein sequence alignment (level 1) and confined alignment (level 2) on functional information obtained from level1 and also potentially from experimental sources.
4.1. Reliability of the data
Multiple resources for obtaining a consensus stretch of the length of the TMDs are used [50]. Sequence alignment programs such as BLAST are widely used to align identified proteins as well as nucleotides with each other to obtain functional information about the query sequence and with it the protein. BLAST shows good performance by combining speed of processing with reliability of the data [51]. The use of one sequence alignment software is compensated by the use of different lengths of amino acid sequences to probe the stability of the data.
4.2. Sequence – structure relation
Predicting the stretch of the TMDs, a large number of programs are used. These programs vary in the algorithm applied to propose the stretches and thus provide proposals which might differs considerably from each other. A threshold of more than 50 % of the programs propose the TMD stretch, is chosen. In addition, a criterion is to allow for a stretch to be long enough to span the bilayer. In some cases, the alignment overlap of all the programs is so that the p-TMDs are almost similar to the r-TMDs predicted by others. In one case the overlap is rather poor, a finding which is also observed for the TMDs of other membrane proteins [50]. It has been argued that this suggest anomalous behavior in terms of the folding of the stretch and consequently on the mechanics and function of the protein. In this respect, TMD1 of NS4B is to some respect close to a region where it is reported that the protein could flip through the membrane [15],[17],[52].
4.3. Sequence–function–relation
The relation is based on the correlation between amino acid sequence and the biological function followed by structural conservation, or phrased differently, that structural evolution is slower than the amino acid evolution [53]. The same functional and structural feature can be anticipated, if the amino acid sequence most similar. Still, non-homologous proteins in respect to their amino acid sequence can adopt the same structure as long as a particular structural motif is needed to perform a specific task [54]. This is particularly the case for the myoglobin / hemoglobin relation where both proteins share very little sequence identity but a lot of structural similarity [55]. It is anticipated that this is in order to take care of the potential discrepancy between amino acids and structural features by using the protocol of using sequence alignment to identify common features amongst proteins for which structural or functional data are available. At this end, the bioinformatic approach identifies a correlation of the NS4B sequence with sequences of proteins of which their functions are found to be involved in regulation of Ca-homeostasis. Thus, the presented data support experimental findings [20].
In addition, similarity of the TMDs 1 and 3 of NS4B with the TMDs of Nsp3 proteins proposes a similar function for NS4B of being involved in the formation of the replication vesicles. Combining the findings, involvement in Ca-homeostasis and vesical formation, Ca-regulation seems to be important for the latter.
5.
Conclusions
Based on these matching results, the proposal of TMDs 2 and 4 playing an active part in homeostasis either as being part of a channeling structural feature or being part of a transporting mechanism, seems plausible. In addition, parts of the stretch containing TMD1 and TMD3 are potentially involved in the formation of replication vesicles.
Experimentally the involvement of NS4B into Ca-homeostasis is proposed, and it would also have been proposed by using the bioinformatics approach. Structural prediction shows the ambiguity for the first TMD, and this sequence stretch in particular, as the programs vary quite substantially in predicting a TMD stretch. Overall, a bioinformatics approach gives thorough evidence of the function and structure of a membrane concerning its transmembrane region.
Based on secondary structure prediction and a two-level sequence alignment approach, the functional role of the TMDs of NS4B is proposed. This algorithm is readily transferable to other proteins and when automatized a valuable tool for predicting structure – function correlations.
DownLoad:
