Elucidation of the interactions between SARS-CoV-2 Spike protein and wild and mutant types of IFITM proteins by in silico methods

COVID-19 is a viral disease that has been a threat to the whole world since 2019. Although effective vaccines against the disease have been developed, there are still points to be clarified about the mechanism of SARS-CoV-2, which is the causative agent of COVID-19. In this study, we determined the binding energies and the bond types of complexes formed by open (6VYB) and closed (6VXX) forms of the Spike protein of SARS-CoV-2 and wild and mutant forms of IFITM1, IFITM2, and IFITM3 proteins using the molecular docking approach. First, all missense SNPs were found in the NCBI Single Nucleotide Polymorphism database (dbSNP) for IFITM1, IFITM2, and IFITM3 and analyzed with SIFT, PROVEAN, PolyPhen-2, SNAP2, Mutation Assessor, and PANTHER cSNP web-based tools to determine their pathogenicity. When at least four of these analysis tools showed that the SNP had a pathogenic effect on the protein product, this SNP was saved for further analysis. Delta delta G (DDG) and protein stability analysis for amino acid changes were performed in the web-based tools I-Mutant, MUpro, and SAAFEC-SEQ. The structural effect of amino acid change on the protein product was made using the HOPE web-based tool. HawkDock server was used for molecular docking and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) analysis and binding energies of all complexes were calculated. BIOVIA Discovery Studio program was utilized to visualize the complexes. Hydrogen bonds, salt bridges, and non-bonded contacts between Spike and IFITM protein chains in the complexes were detected with the PDBsum web-based tool. The best binding energy among the 6VYB-IFITM wild protein complexes belong to 6VYB-IFITM1 (-46.16 kcal/mol). Likewise, among the 6VXX-IFITM wild protein complexes, the most negative binding energy belongs to 6VXX-IFITM1 (-52.42 kcal/mol). An interesting result found in the study is the presence of hydrogen bonds between the cytoplasmic domain of the IFITM1 wild protein and the S2 domain of 6VYB. Among the Spike-IFITM mutant protein complexes, the best binding energy belongs to the 6VXX-IFITM2 N63S complex (-50.77 kcal/mol) and the worst binding energy belongs to the 6VXX-IFITM3 S50T complex (4.86 kcal/mol). The study suggests that IFITM1 protein may act as a receptor for SARS-CoV-2 Spike protein. Assays must be advanced from in silico to in vitro for the determination of the receptor-ligand interactions between IFITM proteins and SARS-CoV-2.


INTRODUCTION
Coronaviruses are crown-shaped, single strand, positive polarity, and enveloped RNA viruses. The virus that causes COVID-19 is classified as beta coronavirus and belongs to the same subgroup as SARS. Therefore it is called SARS-CoV-2 by International Virus Taxonomy Committee [1][2][3]. SARS-CoV-2 interacts with angiotensin converter receptor-2 (ACE2) and binds to the cell. It uses Spike protein to binds the human lung epithelial cells. Spike protein has three chains that are A, B, and C; respectively. Spike protein plays a key role in SARS-CoV-2 infection in human cells with these three chains. It is a major antigene for attachment of the virus via ACE2 receptor to the cell [4].
Spike protein consists of two subunits. S1 subunit is responsible for host cell receptor binding and the S2 subunit is responsible for the fusion of host cell membranes. S1 domain binds the receptors and it functions as the receptor-binding domain (RBD) [5,6]. Spike protein RBD provides the Coronaviruses binding to the host human receptors. Two types of Spike protein have been identified. The structures of the SARS-CoV-2 Spike protein were determined by cryogenic electron microscopy and it finds out two different (open and closed) conformations of the receptor-binding site. The open form is called 6VYB and the closed-form is called 6VXX. The SARS-CoV-2 Spike protein must be in the open conformation to interact with ACE2 to initiate viral entry. This open type triggers severe infection [7,8].
Polymorphisms are genetic differences that occur with a frequency greater than 1% in a population [9]. Polymorphisms are not the cause of diseases but may be a cause of susceptibility to disease. The most common type of polymorphism in the human genome is single nucleotide polymorphisms. SNP is the difference of a single nucleotide in the DNA sequence. SNPs have a major effect on disease susceptibility and drug response. If there is a mutation such as SNPs, structural variations, chromosome number, and structural anomalies in an individual; these mutations or anomalies will change protein structure and function [10]. As a result, many diseases occur. In addition, changes in protein structure also affect the response to drugs [11]. Single nucleotide polymorphisms occurring in the structure of the Spike protein cause changes in the structure of the protein. In this case, small changes in conformation and sequence can lead to differences in the affinity of Spike proteins for ACE2 protein [12].
Interferon-induced transmembrane (IFITM) proteins are a family of small membrane proteins with immunological and developmental roles. Five types of IFITM proteins have been identified in humans: IFITM1, IFITM2, IFITM3, IFITM5, and IFITM10 [13]. The interferon-derived transmembrane protein family (IFITM/MIL/Fragilis) genes encode cell-surface proteins that can regulate cell adhesion and influence cell differentiation. IFITM1 and IFITM3 proteins expressed in primordial germ cells (PGCs) are involved in germ cell development, but their specific functions are unclear [14]. IFITM1 and IFITM3 proteins expressed in PGC are involved in germ cell development, but their specific functions are unclear [ IFITM proteins are a cofactor for SARS-CoV-2 infection and these are potential targets for therapeutic approaches. Peptides and targeting antibodies that derived from IFITM (especially IFITM1, 2, and 3) have important roles in inhibition of SARS-CoV-2 entry and replication in the few cell types as the host human lung cells, cardiomyocytes, and gut [20]. In this study, we demonstrated the interactions (binding energies of the protein-protein complexes and bond types) of wild and mutant types of IFITM1, IFITM2, and IFITM3 proteins with open and closed forms of SARS-CoV-2 Spike protein by molecular docking approach. The workflow chart of the study can be seen in Figure 1. tools were used for the first analysis to determine the effects of the missense SNPs. SIFT is a program that predicts the functional effect of an amino acid substitution in a particular protein. If the score provided by the program is <= 0.05, the effect is deleterious; the exact opposite, the score is > 0.05 then the effect is determined as tolerated [29]. The rsIDs of missense SNPs extracted from NCBI dbSNP were pasted into the SIFT4G predictions tool. Besides the effects of amino acid substitutions, SIFT also provides the information of nucleotide changes and amino acid changes as a table caused by missense SNPs. PROVEAN is a webbased tool in which functional effects are determined based on the amino acid sequence of a protein found in any organism. The default score threshold of the PROVEAN Protein tool is -2.5. In the binary classification, if the score is <= -2.5, the effect is deleterious; the score is > -2.5, the effect is neutral [30]. Amino acid FASTA sequences of IFITM1, IFITM2, and IFITM3 proteins were copied from UniProt (https://www.uniprot.org/) [31] and inserted into the PROVEAN Protein tool with substitutions. PolyPhen-2 works similarly to PROVEAN. PolyPhen-2 scores differ from 0.0 (benign) to 1.0 (damaging) [32]. SNAP2 is a classifier based on a neural network and utilizes evolutionary information from multiple sequence alignments, predicted structural properties, etc. of the amino acid sequences [33]. Mutation Assessor and PANTHER cSNP tools estimate the effects of the substitutions in the evolutionary base. In the first analysis, if the prediction is possibly damaging (POD), probably damaging (PRD), deleterious (DLT), and/or the functional impacts are low, medium, or high in at least 4 of 6 tools, they were saved for the further analysis. Neutral or benign effects were not considered in the study (Table 1).

Data extraction and evaluation of the missense
DDG determination and structural Analysis of the pathogenic forms of IFITM1, IFITM2, and IFITM3 proteins DDG (or referred to as ΔΔG, double changes in Gibbs free energy)s and stability changes of the proteins were determined in the web-based prediction tools: I-Mutant (http://gpcr2.biocomp.unibo.it/cgi/predictors/I-Mutant3.0/I-Mutant3.0.cgi) [34], MUpro (http://mupro.proteomics.ics.uci.edu/) [35] and SAAFEC-SEQ (http://compbio.clemson.edu/SAAFEC-SEQ/#started) [36]. Default settings (Temperature: 25 0 C, pH: 7.0) and the "DDG Value and Binary Classification" option were preferred in I-Mutant. Sequencebased analyses were performed for all pathogenic proteins in the aforementioned tools. If two or three of I-Mutant, MUpro, and SAAFEC-SEQ give the same result, the mutation is considered to make the protein more stable or unstable (Table 2). Also, the HOPE tool was used for structural analysis in terms of size, charge, and physical state of the proteins [37] (Table 3).

Homology modeling of wild and mutant Types of IFITM1, IFITM2, and IFITM3 proteins
Amino acid sequences (FASTA format) of wild-type protein products of IFITM1, IFITM2, and IFITM3 were obtained from UniProt [31]. UniProt IDs for IFITM1, IFITM2, and IFITM3 proteins are P13164, Q01629, and Q01628, respectively. Homology modeling of these proteins was performed via the I-TASSER online server. I-TASSER provides the Confidence (C)-score, estimated TM-score, and estimated RMSD value for each model it creates. The C-score takes a value in the range of [-5,2] and a high C-score indicates high confidence for the related model [38,39]. Also, SWISS-MODEL QMEAN (Qualitative Model Energy Analysis) [40] and ERRAT [41] scores were used for structure assessment. The quality of the model is low when the QMEAN score is less than -4 [42]. If the ERRAT score is higher than 50, the model can be accepted to be of high quality [43]. Mutant-type proteins were generated with the UCSF Chimera program [44]. Simply the wild-type residues were replaced with the mutant-type residues by entering a specific command (swapaa). Swapaa utilizes the information from a rotamer library for this operation [45].

Molecular docking of protein-protein complexes
Molecular docking was carried out between the open and closed forms of Spike proteins of SARS-CoV-2 and wild and pathogenic mutants of the IFITM1, IFITM2, and IFITM3 proteins. For this purpose, the HawkDock server was used for the prediction and analysis of protein-protein complexes [46]. The PDB ID and chain information of the open and closed forms of the Spike proteins were entered into the server (6VYB:A for open form, 6VXX:A for closed-form) as the receptors. PDB files of wild and mutant types of IFITM1, IFITM2, and IFITM3 proteins were uploaded to the server as ligands. Hawkdock server proposes that the larger proteins be designated as receptors (http://cadd.zju.edu.cn/hawkdock/). After molecular docking is complete, Hawkdock presents 100 models for each complex and the scores of the ten of these models with the most negative scores. At this stage, the PDB files of the complexes with the most negative score were saved. MM/GBSA analyzes were performed in the same software and the binding energies were calculated for each complex. The interactions in the protein-protein complexes were visualized with the BIOVIA Discovery Studio 2021 program [47] (Figure 2, Figure 3, and Figure 4) and the H-bond option was selected for the surface display. In addition, all bonds (H-bonds, salt bridges, and non-contact bonds) and the number of these bonds in protein-protein complexes were revealed with the PDBsum online tool (https://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html) [48].

Data extraction and evaluation of the missense SNPs of IFITM1, IFITM2, and IFITM3
There were 1426 SNPs for IFITM1, 1834 SNPs for IFITM2, and 1584 SNPs for IFITM3 in NCBI dbSNP on August 2021. Total missense SNP numbers for IFITM1, IFITM2, and IFITM3 were 95, 139, and 130 respectively. rsIDs of the SNPs, allele changes, amino acid changes, pathogenicity/functional impact scores, predictions, functional impacts, and preservation times according to information from SIFT, PROVEAN, PolyPhen-2, SNAP2, Mutation Assessor, and PANTHER cSNP were shown in Table S1. 18 SNPs have harmful effects on the protein products of IFITM1, IFITM2, and IFITM3 in at least four of the six tools. Among them, the most preserved are rs11552366 (IFITM1), rs200528039 (IFITM1), and rs371179996 (IFITM2) (456 million years). If the preservation is strong, the functional impact of the protein resulting from the variation will be high. Based on the information from whole tools, the most pathogenic missense SNP appears to be rs11552366 (IFITM1) ( Table 1).

DDG determination and structural analysis of the pathogenic forms of IFITM1, IFITM2, and IFITM3 proteins
According to the evaluation results of I-Mutant, MUpro, and SAAFEC-SEQ; amino acid changes in only two (H57L and H36L) of the 18 mutant proteins were predicted to increase the stability of the proteins ( Table 2). Positive DDG values cause mutation to make the protein more stable, while negative values lead to more unstable proteins. The HOPE results showed that the mutant residue in seven (I87T, L41P, P34S, N63S, H3Q, H57L, and H36L) of the 18 mutant proteins had a smaller size compared to the wild-type residue; the mutant residues of the other 11 mutant proteins (Q120H, G74R, E30K, H6Y, V61M, V24M, C84Y, V53M, E20K, S50T, and V90I) were bigger than the wildtype residues. HOPE (Table 3) specified that changes in residue size and hydrophobicity due to mutant residue may affect contact with the lipid membrane. A mutant residue bigger than the wild-type residue may cause bumps. The results indicate that H6Y, N63S, H57L, H36L cause more hydrophobic protein structure. This is in agreement with the I-Mutant and MUpro results, as hydrophobicity can impart stability to mutant proteins. I87T, G74R, C84Y, and P34S lead to the formation of less hydrophobic protein structures, which can make the protein less stable. The reduced hydrophobicity of the protein means loss of H-bonds, resulting in loss of hydrophobic interactions, it also has a negative effect on protein folding. These results are also compatible with the protein structure analysis tools. Table x also shows that the charges of some mutant residues are different from the charges of wildtype residues. The difference in the charge of the mutant-type residue from that of the wild-type residue can cause repulsions within the protein itself or between the protein and the ligand.

Homology modeling of wild and mutant types of IFITM1, IFITM2, and IFITM3 proteins
As a result of homology modeling evaluation for wild-type protein models provided by I-TASSER, the models with the highest C, ERRAT, and QMEAN scores were selected, therefore the homology models with the most quality are Model 1 for IFITM1, Model 1 for IFITM2, and Model 2 for IFITM3 ( Table 4). The C and ERRAT scores appear to be within acceptable ranges, but the QMEAN scores are lower than expected in all models.
Bond analysis (H-bond, salt bridge, and non-contact bond) by PDBSum (Table 6 and Table 7) shows that the number of H-bonds between 6VYB and IFITM1/IFITM3/IFITM3 (wild) is equal (two Hbonds). Hydrogen bonds play an important role in determining the three-dimensional structures and the specificity of ligand binding [49,50]. It is reported that H-bonds displace protein-bound water molecules and promote the ligand-binding affinity [51]. The 6VYB-IFITM3 (wild) complex differs from other 6VYB-IFITM (wild) complexes in that it has a salt bridge in between. In terms of nonbonded contact numbers, there are 158 bonds between 6VYB and IFITM2 (wild), while 150 bonds with IFITM3 (wild) and 131 bonds with IFITM3 (wild). Although the number of non-bonded contacts is less in the 6VYB-IFITM3 (wild) complex than other 6VYB-IFITM (wild) complexes, the extra salt bridge may explain the binding energy of the complex being better than the binding energy of the 6VYB-IFITM2 (wild) complex. The difference between the binding energies of the 6VXX-IFITM (wild) complexes emphasized the importance of a large number of non-bonded contacts. There are 197, 144 and 147 non-bonded contacts between protein-protein in the 6VXX-IFITM1 (wild), 6VXX-IFITM2 (wild) and 6VXX-IFITM3 (wild) complexes, respectively. Although the 6VXX-IFITM3 complex has six H-bonds and four salt bridges, the binding energy of the complex is less negative than the 6VXX-IFITM1 (wild) complex. Non-covalent interactions are relatively weak interactions, but their excess may provide a better interaction than H-bonds and salt bridges.
In Spike-IFITM (wild) complexes, the places where H-bonds are formed are different according to the closed and open forms of the Spike. H-bond formation with the S1 domain of the spike is seen in all 6VYB-IFITM (wild) complexes. Only 6VYB-IFITM1 (wild) complex also has H-bond formation of IFITM1 (wild) with S2. There is H-bond formation with both the S1 and S2 in the 6VXX-IFITM3 (wild) complex, while there is only H-bond formation with the S1 in the 6VXX-IFITM1 (wild) complex. There is no H-bond between 6VXX and IFITM2 (wild). Interesting results have emerged regarding the sites where H-bonds are formed, both with the 6VYB-IFITM (wild) and 6VXX-IFITM (wild) complexes. In the 6VYB-IFITM1 complex, binding occurred between the S1 and extracellular region of IFITM1 and between S2 and the cytoplasmic region of IFITM1. This may indicate that IFITM1 protein acts as a receptor for the Spike protein.
The most negative binding energy belonged to the 6VYB-IFITM1 P34S between the mutant type proteins of IFITM1 and the 6VYB complexes; and showed a little better binding than the wild type of IFITM1. This complex contains two H-bonds at the same positions as the 6VYB-IFITM1 (wild) complex and does not contain salt bridges. This complex contains H-bonds at the same positions as the 6VYB-IFITM1 (wild) complex and does not contain salt bridges. The binding between amino acid 850 (GLU) of S2 and amino acid 36 (HIS) of IFITM1 in one of the cytoplasmic domains in the region may be somewhat restricted by proline, amino acid 34 of wild-type IFITM. Because HOPE results indicated that the new residue (serine) was smaller at this location. In this case, it can be clarified why mutant-type binding is better. There is no 6VXX-IFITM1 (mutant) complex that has a better binding energy value than the 6VXX-IFITM1 (wild) complex.
The number of non-bonded contacts in all of the complexes formed by the mutant type protein of 6VYB and IFITM2 is higher than the non-bonded contacts in the 6VYB-IFITM2 (wild) complex, and therefore the binding energies of all complexes are slightly more negative than the binding energy of the 6VYB-IFITM2 (wild) complex. In the complexes formed by mutant-type proteins of IFITM2 with both 6VYB and 6VXX, the most negative binding energy was observed in the complexes with IFITM2 N63S. The binding energies of 6VYB-IFITM2 N63S and 6VXX-IFITM2 N63S complexes are -25.37 kcal/mol and -50 kcal/mol, respectively.
Among the complexes formed by 6VYB and IFITM3 mutants, the most negative binding energy belongs to the 6VYB-H57L complex (-44.09 kcal/mol). There may be two reasons for this. Compared to the 6VYB-IFITM3 (wild) complex, the number of H-bonds in the 6VYB-H57L complex is 2 times higher (4) and the number of non-bonded contacts is slightly more than double (296). Between the complexes formed by 6VXX and IFITM3 mutants, there is no complex with more negative binding energy than the 6VXX-IFITM3 (wild) complex. Between the complexes formed by 6VXX and IFITM3 mutants, there is no complex with more negative binding energy than the 6VXX-IFITM3 (wild) complex. The 6VXX-IFITM3 (wild) complex has more H-bonds and salt bridges than the 6VXX-IFITM3 mutant complex. An interesting result here is that the binding energy (4.86 kcal/mol) of the 6VXX-S50T complex is positive, making it the worst binding. Unlike other 6VXX-IFITM3 complexes, in the 6VXX-S50T complex, three H-bonds are formed with the first and second amino acids of the signal peptide located at the N-terminus of the Spike protein and the cytoplasmic domain of IFITM3 S50T. It is difficult to say whether this is related to any mechanistic event, since there may also be a binding affinity seen only in silico.

CONCLUSION
Although effective vaccines have been developed for COVID-19, there are still points to be clarified about the mechanism of SARS-CoV-2. The mechanisms need to be well known in the development of new inhibitory drugs for the treatment of viral diseases. This study is the only study that offers a molecular docking approach to determine the interaction of both wild and mutant types of IFITM proteins with the Spike protein of SARS-CoV-2. With this study, we propose for the first time that there may be another receptor (IFITM1) for Spike other than ACE2, TMPRSS2 and FURIN. Our study is in silico-based and needs to be investigated in vitro as well.  [21] Home -SNP -NCBI, (n.d.