Bioinformatics Analysis of Domain 1 of HCV-Core Protein: Iran.

Hepatitis C virus (HCV) infection is a serious global health problem and a cause of chronic hepatitis, liver cirrhosis, and hepatocellular carcinoma (HCC). Bioinformatics software has been an effective tool to study the HCV genome as well as core domains. Our research was based on employing several bioinformatics software applications to find important mutations in domain 1 of core protein in Iranian HCV infected samples from 2006 to 2017, and an investigation of general properties, B-cell and T-cell epitopes, modification sites, and structure of domain 1. Domain 1 sequences of 188 HCV samples isolated from 2006 to 2017, Iran, were retrieved from NCBI gene bank. Using several tools, all sequences were analyzed for determination of mutations, physicochemical analysis, B-cell epitopes prediction, T-cell and CTL epitopes prediction, post modification, secondary and tertiary structure prediction. Our analysis determined several mutations in some special positions (70, 90, 91, and 110) that are associated with HCC and hepatocarcinogenesis, efficacy of triple therapy and sustained virological response, and interaction between core and CCR6. Several B-cell, T-cell, and CTL epitopes were recognized. Secondary and tertiary structures were mapped fordomain1 and core proteins. Our study, as a first report, offered inclusive data about frequent mutation in HCV-core gene domain 1 in Iranian sequences that can provide helpful analysis on structure and function of domain 1 of the core gene. © Springer Nature B.V. 2019.

Introduction

n class="Disease">HCV infection is a n class="Chemical">serious global health problem and causes n class="Disease">chronic hepatitis, liver cirrhosis, and HCC (Akuta et al. 2007; Ajorloo et al. 2015; Alborzi et al. 2017, 2015; Moayedi et al. 2018; Hashempoor et al. 2018). It is estimated that 160 million people are infected worldwide (Lavanchy 2011). The prevalence rate of HCV infection is from 0.2 up to 40% in different countries, and this prevalence in Iran is 0.16% (Sefidi et al. 2013).

Lack of an effective vaccipan class="Chemical">ne and therapeutic choices has leaded to the rapid growth of n class="Disease">HCV infection (Lauer and Walker 2001).

n class="Species">HCV has six ln class="Chemical">arge different genotypes (1–6) and 70 distinct subtypes (a, b, c, etc.) globally (Martro et al. 2011). Genotyping analysis showed 30–33% difference in each genotype, and in subtypes around 20–25% (Sefidi et al. 2013). According to the current investigations in Iran, the predominant n class="Species">HCV subtype is 1a, followed by 3a and 1b (Sefidi et al. 2013).

n class="Species">HCV is a positive-strand Rpan> class="Chemical">NA virus encoding three structural components, the core protein, and two E1 and E2 envelope glycoproteins (Ajorloo et al. 2015).

The core protein has many confirmed roles: core binds RNA anpan>d Dpan> class="Chemical">NA and has an important function in RNA packing. It has been determined that HCV-core is a nucleic acid chaperone similar to retroviral nucleocapsid (NC) proteins in act acting to rearrange HCV 3′UTR, resulting in RNA dimerization in vitro (Caval et al. 2011; Cristofari Gl; Ivanyi-Nagy et al. 2004; Steinmann et al. 2008).

n class="Species">HCV-core is a highly basic protein that forms the viral pan> class="Chemical">NC and has interactions with cellular proteins and signal transduction pathways. As a result of HCV-core and host cell interactions, core may have a function in persistent n class="Disease">infection and the pathogenesis of HCV liver disease (Polyak et al. 2006; Hashempour et al. 2015).

Core protein copan class="Chemical">nsisted of n class="Chemical">three predicted domains: aa1-aa117 domainclass="Chemical">pan> 1 (domain D1), aa117-aa177 domain 2 (domain D2), n class="Gene">and 177–191 domain 3 (domain D3) (Strosberg et al. 2010).

Domain 1 contains fren class="Chemical">quenpan>t positively chpan> class="Chemical">arged amino acids, and is involved in RNA binding, promotes dimerization of the viral RNA, and has a significant role in NC formation and core envelopment by endosomal membranes (Ivanyi-Nagy et al. 2006).

Several identified mutations in domain 1 are involved in the development of n class="Disease">HCC and hepatocarcinogenesis, the efficacy of triple therapy, anpan>d interactionpan> betweenpan> core anpan>d pan> class="Gene">CXCL6 (Akuta et al. 2007, 2010, 2011; Fishman et al. 2009; Ogata et al. 2002; Takahashi et al. 2001; Idrees and Ashfaq 2013).

Humoral and cellular immune responses against n class="Disease">HCV infections are inefficienpan>t anpan>d there is no conpan>vincing explanpan>ationpan> to unpan>derstanpan>d pan> class="Species">HCV immune pathogenesis (Gremion and Cerny 2005). However, HCV-specific IgM and IgG together were detected in acute infection. There are some outstanding proofs supporting a role for Abs in control of HCV infection and especially in reinfection (Cashman et al. 2014).

Some researchers have described a rise in anti-n class="Species">HCV humeral immunpan>e after immunclass="Chemical">pan>ization with core protein (Aghasadeghi et al. 2006). Other researchers have claimed that in pan> class="Disease">acute HCV infection the titer of antibodies is very low, and delay in neutralizing antibody production is responsible for ineffective ability to prevent HCV infection (Netski et al. 2005). Other sources have introduced a large diversity of epitopes in n class="Species">HCV proteins as a possible way to escape the humoral response (Pavio and Lai 2003).

Cellular immune responses have an important role in the clearance of n class="Disease">HCV infection. pan> class="Species">Patients with cellular immune dysfunction like human immunodeficiency virus (HIV) have rapid HCV progression. Some researchers have claimed that this system promotes liver injury by cytolysis activity of infected cells. In spite of cellular immune responses, HCV often evades recognition and has the ability to persist (Ward et al. 2002).

Several studies have described a number of pathways for n class="Species">HCV to escape from cell responclass="Chemical">pan>ses, including impaired oligo-/monpan>o-specific or no virus-specific pan> class="Gene">CD4+ and CD8+, mutation of epitopes, weakness of proliferative capacity, and cytotoxicity and ability to secrete TNF-α and IFN-γ by n class="Gene">CD8+ T cells (Neumann-Haefelin et al. 2005). In addition, regulatory T cells (Tregs) induced by HCV infection have a significant role in the impaired activity of cellular immune response (Hashempour et al. 2015; Hashempoor et al. 2010).

Bioinformatics tools are efficient means to study viruses and different parts of the n class="Species">HCV genpan>ome like core domains (Idrees anpan>d Ashfapan> class="Chemical">q 2013; Moattari et al. 2015; Dehghani et al. 2017; Nezafat et al. 2018; Atapour et al. 2018; Sarvari et al. 2014; Behzad Dehghani and Zahra Hasanshahi 2019). Bioinformatics tools are efficient means to study viruses and different parts of the HCV genome like core domains. Many programs have been developed to analyze function, structures, and modification of core protein, providing a large amount of information about core domains and important mutation sites (Akuta et al. 2007, 2010, 2011).

Current data are useful for predictiopan class="Chemical">n of n class="Disease">HCV disease development anclass="Chemical">pan>d treatment responpan>se. Inpan> this study, we employed several bioinformatics tools to find important mutationpan>s in domain 1 of the core protein, general properties of B-cell and T-cell epitopes, modificationpan> sites, and structure of domain 1 in n class="Disease">Iranian HCV infected samples from 2006 to 2017.

Materials and Methods

Sequence Alignment and Phylogenetic Tree

Domain 1 sen class="Chemical">quenclass="Chemical">pan>ces of 188 Iranianpan> pan> class="Species">HCV samples and reference sequences of HCV genotypes that were registered in n class="Chemical">NCBI gene bank (http://www.ncbi.nlm.nih.gov/) from 2006 to 2017 were downloaded. Homology among sequences was determined using multiple sequence alignment available in CLC- sequence viewer software under the following parameters: gap open cost, 10; gap extension cost, 1.0; and very accurate progressive alignment algorithm. Also, phylogenetic trees were analyzed through CLUSTAL X software, version 1.81, by neighbor-joining times to confirm the reliability of phylogenetic trees. The accession numbers of all sequences are displayed in Table 1.

Table 1

The accession numbers of all 188 sequences that were used in this study

2006(16)	DQ065821-DQ065836 and DQ202323
2012(1)	JN129986
2013(84)	KC118250-KC118333, KC118288
2014(28)	KC285337-KC285362, JQ341409,KF218585,KC118320
2016(59)	KT329291-KT329341, KU736828-KU736837

The accession pan class="Chemical">numbers of all 188 sen class="Chemical">quences that were used inclass="Chemical">pan> this study

Determination of Mutations

By considering previous studies, several significant mutations that are involved in: n class="Disease">1-hepatocellular carcinoma, 2-viral responpan>se to triple therapy, anpan>d 3-the interactionpan> betweenpan> core anpan>d class="Chemical">pan> class="Gene">CXCL6, were determined. All sequences were compared to find mentioned mutations (Akuta et al. 2007; Fishman et al. 2009; Ogata et al. 2002).

Physico-chemical Analysis

General properties of domain 1 (genotype 1a) were determined by employing “Expasy’sProtParam” (http://expasy.org/tools/protparam.html) and ProtScale at (http://web.expasy.org/protscale/) (Gasteiger et al. 2005).

B-Cell Epitopes Prediction

Chou and Fasmapan class="Chemical">n, Karplus and Schulz, Kolaskar & Tongaonkar, Emini, Parker, and BepiPred methods at http://www.immuneepitope.org (http://tools.immuneepitope.org/tools/bcell/iedb_input) were run for prediction of B-Cell epitopes positions (Chou and Fasman 2009; Karplus and Schulz 1985; Emini et al. 1985; Parker et al. 1986; Larsen et al. 2006).

On hydrophilicity, flexibility/mobility, accessibility, polarity, exposed surface and turns features by BcePred (http://www.imtech.res.in/raghava/bcepred) B-cell epitopes prediction were performed (Saha and Raghava 2004).

ABCpred software (http://www.imtech.res.in/raghava/abcpred/) predicted 16 meric B-cell epitopes (Saha and Raghava 2006a, b).

Prediction of T-Cell, CTL Epitopes and Allergic Properties

ProPred-I (http://www.imtech.res.in/raghava/propred1/) (Singh and Raghava 2003) was employed for n class="Gene">MHC Class-I binding peptide predictionpan> anpan>d proposed (http://www.imtech.res.inclass="Chemical">pan>/raghava/propred/) was used for n class="Gene">MHC Class-II binding peptide prediction. Programs were worked at a 4% default n class="Chemical">threshold by the proteasome and immunoproteasome filters on at 5% threshold (Singh and Raghava 2001).

n class="Gene">MHC class I and IIpan> predictions were determined using the Immune Epitope Database (IEDB) (http://tools.immuneepitope.org/main/). For predictionpan> of CTL epitopes, “ctlpred” and ANN methods were used (48).

Probability of apan class="Chemical">ntigenicity was expected by VaxiJen software at http://www.ddg-pharmfac.net (Doytchinova and Flower 2007).

IgE epitopes and papan class="Chemical">n class="Disease">allergic properties were estimated at http://www.imtech.res.inclass="Chemical">pan>/raghava/algpred/index.html by using AlgPred (Saha and Raghava 2006).

Post-modification

n class="Chemical">Serine, n class="Chemical">threonine, and n class="Chemical">tyrosine phosphorylation sites prediction was done using DISPHOS (http://www.dabi.temple.edu/disphos/pred.html) (Iakoucheva et al. 2004) and NetPhos (http://www.cbs.dtu.dk/services/NetPhos/) (Blom et al. 1999). Kinase specific phosphorylation sites were determined by NetPhosK (http://www.cbs.dtu.dk/services/NetPhosK/) (Blom et al. 2004). NetNGlyc (http://www.cbs.dtu.dk/services/NetNGlyc/) (Gupta and Brunak 2002) and GlycoEP (http://www.imtech.res.in/raghava/glycoep/submit.html) were employed for N-glycosylation sites prediction (Chauhan et al. 2013).

Secondary and Tertiary Structure Prediction

To predict secondary and tertiary structures of core and domain 1 of genotype 1a, SOPMA at (http://npsa-pbil.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma.html) (Geourjon and Deleage 1995), I-TASn class="Chemical">SER at (http://zhanpan>glab.ccmb.med.umich.edu/I-TASpan> class="Chemical">SER) (Roy et al. 2010), PHYRE2server at (http://www.sbg.bio.ic.ac.uk/~phyre2/html) (Kelley and Sternberg 2009), (PS)2-v2 Server at (http://ps2v2.life.nctu.edu.tw) (Chen et al. 2006) were employed. To evaluate the stereochemistry and quality of 3D structures Qmean at (http://swissmodel.expasy.org/qmean/cgi/index.cgi) (Benkert et al. 2008) was used, and the Ramachandaran plot was mapped by Rampage (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php).

The Signal Peptide Prediction

The Signal peptide was predicted by “Sigpan class="Chemical">nal-BLAST”, and “SignalP 4.1 n class="Chemical">Server”.

Prediction of Epitopes Digestion: Peptide Cutter was Used to Determine Potential Cleavage Sites Cleaved by Proteases

Research Ethics

All data were collected anonymously in accordance with legal ren class="Chemical">quiremenpan>ts regarding data protectionpan> anclass="Chemical">pan>d medical conpan>fidenpan>tiality. Approval from the Faculty pan> class="Species">Human Research Ethics Committee (Shiraz University of Medical Sciences) was obtained before the commencement of the study.

Results

Mutation and Phylogenic Tree

By considering all submitted sen class="Chemical">quenpan>ces inclass="Chemical">pan> n class="Chemical">NCBI GenBank we could not find any sen class="Chemical">quences related to 2017.

Phylogenetic tree for all sen class="Chemical">quenpan>ces was shownclass="Chemical">pan> in Fig. 1. All 2006 sepan> class="Chemical">quences were placed in a cluster at the bottom of the tree, and a sequence of 2012 has a high similarity to KF218585.1 (2014). The majority of sen class="Chemical">quences were closer to 1a and 3a than other reference sequences.

Fig. 1

Phylogenetic tree based on domain1 sequences and by using neighbor joining method. The phylogenetic tree was constructed by the NJ method. The numbers at the forks show the numbers of occurrences of the repetitive groups to the right out of 100 bootstrap samples. All used reference sequences were showed after accession numbers (1a, 1b, and etc.). Sequences were categorized in five major clusters

Phylogenetic tree based on domain1 sen class="Chemical">quenpan>ces anpan>d by using neighbor joining method. The phylogenpan>etic tree was conpan>structed by the pan> class="Chemical">NJ method. The numbers at the forks show the numbers of occurrences of the repetitive groups to the right out of 100 bootstrap samples. All used reference sequences were showed after accession numbers (1a, 1b, and etc.). Sequences were categorized in five major clusters

All important mutatiopan class="Chemical">n positions were listed in Table 2; the majority of mutations happened in 2013 and 2016 samples. n class="Chemical">No mutationpan> was detected in12, 23, 25, 39, 45 positionpan>s.

Table 2

Frequency of all identified mutations in HCV-core domain1 sequences

	2006	2012	2013	2014	2016	Ref
11(T)	–	–	1(P)1.2%	–		1a(N)
12(K)	–	–	–	–		–
23(K)	–	–	–	–		–
25(P)	–	–	–	–		––
39(R)	–	–	–	–		–
45(G)	––	–	–	–		–
49(T)	2(V)11.7%	–	2(P) 2.4%, 5(A) 6.1%, 1(V) 1.2%	1(V) 3.5%	1(p)1.6%	–
69(R)	–	–	–	–		–
70(R)	2(Q)11.7%	–	3(Q) 3.7%, 2(H) 2.4%, 1(P) 1.2%	4(Q)14.3%	25(Q)42.3%,1(H)1.6%	3a,6a (Q)
74(R)	–	–	–	–		–
78(Q)	–	–	–	1 (R) 3.5%		2a,2b (K)
90(G)	–	–	1(L) 1.2%, 1(S)1.2%	–		–
91(C)	2(M)11.7%	–	13(M)16.4%, 3(L) 3.7%, 1(F) 1.2%	–	43(M)72.8%, 7(L)11.8	1b (M); 2a and 5a (L)
102(G)	–	–	–	–		5a (S)
110	8(T)47%, 9(N)53%	–	23(N)28.4%, 5(S)6.1%, 53(T)65.4%	18(T)64.3%, 8(N)28.57%	N (10)17%	1b, 1a, 2b(T); 1c(S); 4a, 3a, 6a, 2a, 5a(N)

The majority of mutations happened in amino acid residues 49, 70, 91, 110

Fren class="Chemical">quency of all idenpan>tified mutationpan>s in pan> class="Species">HCV-core domain1 sequences

The majority of mutations happepan class="Chemical">ned in amino acid residues 49, 70, 91, 110

Protparam analysis

Protparam results for domaipan class="Chemical">n 1 are listed in Table 3. Because of the high percentage of n class="Chemical">basic amino acids, domain 1 is a highly basic peptide (Theoretical pI: 12). The instability index, anclass="Chemical">pan> estimate of the stability of a protein in a test tube, showed that domain 1 is an unstable peptide. Aliphatic index, a positive factor for the increase of thermostability of proteins, indicated that this peptide is a thermostable peptide. GRAVY is a hydropath city index and increasing positive score indicates a greater hydrophobicity, so this peptide is a hydrophilic peptide.

Table 3

Domain1 physicochemical properties computed by “Protparam”

Amino acid composition	Number	Percentage (%)
Number of amino acids: 117
Molecular weight: 13320.2
Theoretical pI: 12.05
Ala (A)	4	3.40
Arg (R)	21	17.90
Asn (N)	5	4.30
Asp (D)	2	1.70
Cys (C)	1	0.90
Gln (Q)	6	5.10
Glu (E)	3	2.60
Gly (G)	16	13.70
Ile (I)	2	1.70
Leu (L)	6	5.10
Lys (K)	7	6.00
Met (M)	1	0.90
Phe (F)	1	0.90
Pro (P)	17	14.50
Ser (S)	7	6.00
Thr (T)	6	5.10
Trp (W)	5	4.30
Tyr (Y)	3	2.60
Val(V)4		3.40
Total number of negatively charged residues (Asp + Glu)	5
Total number of positively charged residues (Arg + Lys)	28
The estimated half-life is
Mammalian reticulocytes, in vitro	30 h
Yeast, in vivo	> 20 h
Escherichia coli, in vivo	> 10 h
The instability index (II)	77.69 unstable
Aliphatic index	40
Grand average of hydropathicity (GRAVY)	− 1.424

Domain1 physicochemical properties computed by “pan class="Chemical">Protparam”

ProtScale Analysis

Hydropathicity analysis by Kyte J. and Doolittle R.F. method showed that the major part of the peptide had a negative score; the maximum hydropath city score was on n class="Chemical">aa34 (valin) anpan>d the minimum hydropath city score was onclass="Chemical">pan> aa 14 (n class="Chemical">asparagine).

Amino acids flexibility predicted by Bhaskarapan class="Chemical">n R Ponnuswamy P.K method indicated that the maximum flexibility was around amino acid 58 (n class="Chemical">proline) and the minimum was around aa 95 (n class="Chemical">glycine).

Transmembrane (TM) tendency calculated by Zhao, G., London E. method, showed that the major part of peptide had a negative score, and the maximum transmembrane tendency was on n class="Chemical">aa34 (pan> class="Chemical">valine) and the minimum was on aa10 (lysine).

Peptide polarity predicted by Grantham R method showed that the maximum polarity was on amino acid 12 n class="Gene">and 13 (pan> class="Chemical">lysine and arginine) and the minimum polarity was on aa 34 (valine).

http://www.immuneepitope.org opan class="Chemical">nline software:

Chou and Fasman Beta-Turn Prediction, which is based on the rationale for predicting turns to predict antibody epitopes, showed one high score region, 99–112.

Emini Surface Accessibility pan class="Chemical">Prediction, which is based on surface accessibility scale, showed two high score positions (4–20 and 49–61).

Karplus and Schulz flexibility scale was used for B-cell prediction; this method is based on mobility of protein segments on the basis of the known temperature B factors of the a-n class="Chemical">carbons of 31 proteins of knpan>ownpan> structure. Results demonpan>strated two positionpan>s with the highest score (50–60 anpan>d 5–12).

Kolaskar & Tongaonkar Antigenicity method is based on physicochemical properties of amino acid residues and their fren class="Chemical">quenpan>cies of occurrenpan>ce in experimenpan>tally knpan>ownpan> segmenclass="Chemical">pan>tal epitopes to predict antigenpan>ic determinanpan>ts onpan> protein. Results showed five positionpan>s of 20–39, 43–49, 63–69, 78–85, anpan>d 93–101.

Parker Hydrophilicity Prediction method is based on peptide retention times during high-performance lin class="Chemical">quid chromatography (HPLC) onpan> a reversed-phase column. By this method, two regionpan>s were founpan>d: 9–16, 51–57.

Linear B-cell epitopes were determined using BepiPred. This method is based on a combination of a hidden Markov model and a propensity scale method. n class="Chemical">Three regionpan>s (1–25, 51–84, pan> class="Gene">and 102–116) were founded by BepiPred analysis.

Bcepred Results

For a combinatiopan class="Chemical">n of all physicochemical properties (hydrophilicity, flexibility/mobility, accessibility, polarity, exposed surface, and turns) for linear B-cell epitope prediction based on physicochemical properties on a non-redundant dataset:

Using bcepred opan class="Chemical">nline software, n class="Chemical">three regionpan>s (4–22, 47–59, n class="Gene">and 109–115) with the highest combined score were found.

Five 16 meric conpapan class="Chemical">n class="Chemical">served regionpan>s (78, 49, 1, 64 n class="Gene">and 102) were found by ABCpred prediction n class="Chemical">Server (Table 4).

Table 4

16 meric conserved B-cell epitopes regions in HCV-core domain1, predicted by ABCpred online software

Rank	Sequence	Position	Score
1	QPGYPWPLYGNEGCGW	78–93	0.92
2	TRKTSERSQPRGRRQP	49–64	0.9
3	MSTNPKPQKKNKRNTN	1–16	0.89
4	PIPKARRPEGRTWAQP	64–79	0.86
4	GSRPSWGPTDPRRRSR	102–117	0.86

16 meric conn class="Chemical">served B-cell epitopes regionclass="Chemical">pan>s in pan> class="Species">HCV-core domain1, predicted by ABCpred online software

VaxiJen Prediction

According to a predefipan class="Chemical">ned cutoff of VaxiJen program, domain 1 was confirmed as a probable antigen (model: virus and n class="Chemical">threshold: 0.4).

IgE Epitopes

The prediction of allergenic proteins by mapping of IgE epitope, SVM, and hybrid methods showed that domain 1 was not an allergen protein.

T Cell-Epitopes

Regarding T-cell responses against n class="Species">HCV, previous researches founpan>d some hosts’ pan> class="Species">human leukocyte antigen (HLA) alleles associations with HCV infection in Iranian patients. We found several epitopes of HLA’s shown in Table 5.

Table 5

HLA predicted epitopes in HCV-core domain1 sequence (genotype 1a) for HLA’s that were determined by previous researches in Iranian patients

HLA types	Epitopes positions
HLA-A1	1_10
HLA-A2	29–44, 77–85
HLA-A3	2–10,29–59, 96–104
HLA-B*3501	6–14, 27–35, 41–49, 57–65,78–91, 99–116
HLA-B*3801	29–37, 77–85, 89–97
HLA-Cw*0401	23–37, 58–65, 78–93
HLA-Cw*0702	27–35, 73–86
HLA-DQA1*05:01	2–16,13–27
HLA-DQA1*0201	NO
HLA-DQB1*0602	2–16,17–31, 102–116
HLA-DQB1*0301	96–104
HLA-G: HLA-G*01:01, HLA-G*01:02	4–16,8–21, 22–33, 66–79, 100–112
HLA-DRB1*0301	96–104
HLA-DRB1*0305	96–104
HLA-DRB1*0309	96–104
HLA-DRB1*0701	No
HLA-DRB1*11	34–42

n class="Gene">HLA predicted epitopes in pan> class="Species">HCV-core domain1 sequence (genotype 1a) for n class="Gene">HLA’s that were determined by previous researches in Iranian patients

Some studies found both papan class="Chemical">n class="Gene">CD4 helper anclass="Chemical">pan>d n class="Gene">CD8 CTL responses against n class="Disease">HCV infection. ctlpred found several epitopes for CTL (Table 6).

Table 6

CTL epitopes in HCV-core domain 1; the high score epitopes are displayed

Peptide Rank	Position	Sequence	Score	Prediction
1	48–56	ATRKTSERS	1	Epitope
2	80–88	GYPWPLYGN	1	Epitope
3	2–10	STNPKPQKK	0.99	Epitope

CTL epitopes in papan class="Chemical">n class="Species">HCV-core domainclass="Chemical">pan> 1; the high score epitopes are displayed

Postmodification

Predictiopan class="Chemical">n of serine, threonine, and tyrosine phosphorylation sites by DISPHOS showed one position (116) in domain1.

By “n class="Chemical">Npan>etPhos” software we found 10 phosphorylationpan> sites (Fig. 2), 6 sites for serine (53, 56, 99, 103, 106, and 116) 3 sites for threonine (15, 49, and 52) and one site for n class="Chemical">tyrosine (86).

Fig. 2

Phosphorylation sites prediction for domain 1 using “NetPhos” online software. Green lines indicate 6 sites for serine, blues lines show 3 sites for threonine, and one purple line shows tyrosine. All sites with scores above the threshold of 0.5 were considered as phosphorylation sites

Phosphorylation sites prediction for domain 1 using “n class="Chemical">NetPhos” onpan>line software. Greenpan> lines inclass="Chemical">pan>dicate 6 sites for n class="Chemical">serine, blues lines show 3 sites for n class="Chemical">threonine, and one purple line shows tyrosine. All sites with scores above the threshold of 0.5 were considered as phosphorylation sites

n class="Chemical">NetPhosK results determined four phosphorylationpan> sites, three threonine amino acids (3, 15, and 49) for protein kinase C and one n class="Chemical">serine (116) for protein kinase A. No glycosylation site was found by NetNGlyc and GlycoEP.

Secondary structure prediction for core and domain1 by using SOMPA software was summarized in Table 7 and Fig. 3.

Table 7

Percentage of secondary structures in core and domain1

	Alpha helix (Hh)	Extended strand (Ee)	Beta turn (Tt)	Random coil (Cc)
Core	39, 20.42%	25, 13.09%	14, 7.33%	113, 59.16%
Domain 1	NO	19, 16.24%	10, 8.55%	88, 75.21%

Fig. 3

Secondary structure prediction using SOMPA. Red region is extended strand, blue is the alpha helix, green is beta turn, and purple is the random coil. The majority of core structure belongs to random coil

SOMPA showed there was no alpha helix structure in domain1 and the major part of it was the random coil. But the combination of (PS)2-v2 and PHYRE2 showed there was an alpha helix structure in 8–15 region (Figs. 4, 5). All programs displayed extended strand in the 29–36 region.

Fig. 4

Secondary structure prediction using PHYRE2. The result of this tool shows that the majority of the core structure (40%) is alpha helix which is indicated with green helix, also the confidence keys of the predicted structure for these regions are high

Fig. 5

Secondary structure prediction using (PS)2-v2. C coil, H helix, and E extended strand. The majority of core structure contains coil structure

Secondary structure predictiopan class="Chemical">n using (PS)2-v2. C coil, H helix, and E extended strand. The majority of core structure contains coil structure

3D structures were determined by all n class="Chemical">three onclass="Chemical">pan>line software but onpan>ly structures predicted by I-TASpan> class="Chemical">SER were reliable. Final structures (Figs. 6, 7) were validated by Qmean. QMEAn class="Chemical">Nscore and Z-score for calculated for core were 0.242 and − 5.43. The scores were not satisfactory but at least provided an overview of the core protein structure. QMEANscore and Z-score for domain 1 were 0.61 and − 1.24 confirming the quality and reliability of the predicted structure.

Fig. 6

3D structure of the core protein using “I-TASSER” program. The selected model had the highest C score and it was qualified by “Qmean”

Fig. 7

3D structure of domain1 using I-TASSER program. The selected model had the highest C score and it was qualified by “Qmean”

3D structure of the core protein usipan class="Chemical">ng “I-TASSER” program. The selected model had the highest C score and it was class="Chemical">pan> class="Chemical">qualified by “Qmean”

3D structure of domain1 usipan class="Chemical">ng I-TASSER program. The selected model had the highest C score and it was class="Chemical">pan> class="Chemical">qualified by “Qmean”

Ramachandran plot was assessed by RAMPAGE, and percentages of the favoured region, allowed region, and outlier region for core were 63.0%, 27.5%, and 9.5% respectively (Fig. 8). RAMPAGE results for domain1 showed 61.7% of residues in favored region and 24.3 in allowed region (Fig. 9). Figure 10 showed the showed T cell and B-cell epitopes on the surface of the core protein.

Fig. 8

Ramachandran plot was used to visualize energetically allowed regions for backbone dihedral angles ψ against ϕ of amino acid residues in modeled protein structure (LCC model) for tertiary structure of core protein by RAMPAGE; the majority of amino acids residues were in favored region (119 amino acids) and allowed region (52 amino acids)

Fig. 9

Ramachandran plot was used to visualize energetically allowed regions for backbone dihedral angles ψ against ϕ of amino acid residues in modeled protein structure (LCC model) for tertiary structure of domain1 by RAMPAGE; majority of amino acids residues were in favored region (71 amino acids) and allowed region (28 amino acids)

Fig. 10

A: the position of the B-Cell epitopes (yellow region) on core tertiary structure and B: the T-cell epitopes (yellow region) on core 3D structure

Ramachandrapan class="Chemical">n plot was used to visualize energetically allowed regions for backbone dihedral angles ψ against ϕ of amino acid residues in modeled protein structure (LCC model) for tertiary structure of core protein by RAMPAGE; the majority of amino acids residues were in favored region (119 amino acids) and allowed region (52 amino acids)

A: the position of the B-Cell epitopes (yellow regiopan class="Chemical">n) on core tertiary structure and B: the T-cell epitopes (yellow region) on core 3D structure

Core Signal Peptide

Both onlipan class="Chemical">ne tools “Signal-BLAST” and “SignalP 4.1 n class="Chemical">Server” were not able to predict any signal peptide for core protein.

Cleavage Sites Prediction

The results of “PeptideCutter” prediction were summarized in Table 8. The prediction was done for all the predicted epitopes. According to the results, the antigenic epitopes that had the lower number of cleavage positions for enzymes were more potential for B cell or T cell epitopes.

Table 8

The results of predicted cleavage positions for 12 common proteases: B-cell, T-cell, and CTL predicted epitopes

Positions	Caspase	Chymotrypsin	Clostripain	Elastase	Pepsin	Proteinase K	Asp-N	Staphylococcal peptidase I	Thrombin	Proline	Trypsin
B-cell epitopes
78–93		3(8,9,16)			5(4, 7, 9,15, 16)	6(4, 6, 8, 9, 12, 16)	1(11)	1(12)
49–64			5 (2, 7, 11, 13, 14)			3(1, 4, 6)	1(5)	1(6)			6(2, 3, 7, 11, 13, 14)
1_16		1(1)	1(13)			2(3,15)				1(7)	4(9, 10, 12, 13)
64–79		1(13)	3(6, 7, 11)	2(5,14)	1(12)	6(2, 5, 9, 12, 13, 14)	1(8)	1(9)		1(8)	3(4, 6, 11)
102–117		1(6)	5(3, 12, 13, 14, 16)			2(6,9)	1(9)			1(4)	3(12, 14, 16)
CTL epitopes
48–56			2(3,8)	1(1)		4(1, 2, 5, 7)	1(6)	1(7)			3(3, 4 ,8)
80–88		2(6,7)			3(2, 5, 7)	4(2, 4, 6, 7)
2_10						1(2)				1(6)	2(8,9)
T-cell epitopes
HLA-A1
1_10		1(1)				1(3)				1(7)	2(9,10)
HLA-A2
29–44		3(6, 7, 15)	3(10, 11, 14)	2(2, 5)	4(5, 6, 8, 15)	7(1, 2, 5, 6, 7, 8, 15)			1(10)		3(10, 11, 14)
HLA-A2
77–85		1(9)		1(1)	2(5,8)	4(1, 5, 7, 9)
HLA-A3
2_10						1(2)				1(6)	2(8,9)
HLA-A3
29–59	1(25)	3(7,8,16)	7(11, 12, 15, 19, 22, 27, 31)	4 (3, 6, 18, 20)	4(6 7 9 16)	12(2, 3, 6, 7, 8, 9, 16, 18, 20, 21, 24, 26,)		1(26)	1(11)		8(11 12 15 19 22 23 27 31)
HLA-A3
96–104
HLA-B*3501
6_14			1(8)							1(2)	4(4, 5, 7, 8)
HLA-B*3501
27_35		1(9)		2(5,8)	2(8,9)	4(4, 5, 8, 9)
HLA-B*3501
41–49		1(4)	2(3,7)	2(6,8)	1(4)	4(4 6 8 9)					2(3,7)
HLA-B*3501
57–65			3(3 5 6)			1(9)					4(3, 5, 6, 11)
HLA-B*3501
78–91		2(8,9)			3(4, 7, 9)	5(4, 6, 8, 9, 12)	1(11)	1(12)
HLA-B*3501
99–116		1(9)	5(3, 6, 15, 16, 17)			2(9, 12)	1(12)			1(7)	3(3, 15, 17)
HLA-B*3801
29–37		3(7, 8, 9)	2(3, 6)	4(6, 7, 8, 9)		6(2 3 6 7 8 9)
HLA-B*3801
77–85		1(9)		1(1)	2(5,8)	4(1, 5, 7, 9)
HLA-B*3801
89–97		3(5, 8, 9)		1(6)	5(4, 5, 7, 8, 9)	5(1, 5, 6, 8, 9)		1(1)
HLA-Cw*0401
23–37		3(7, 8, 9)		2(3, 6)	4(6, 7, 8, 9)	6(2, 3, 6, 7, 8, 9)
HLA-Cw*0401
58–65			3(2,4,5)								3(2,4,5)
HLA-Cw*0401
78–93		3(8, 9 ,16)			5(4, 7, 9, 15, 16)	6(4, 6, 8, 9, 12, 16)	1(11)	1(12)
HLA-Cw*0702
27–35		1(9)		2(5, 8)	2(8,9)	4(4, 5, 8, 9)
HLA-Cw*0702
73–86		3(4, 13, 14)	1(2)	1(5)	4(3 9 12 14)	7(3 4 5 9 11 13 14)					1(2)
HLA-DQA1*05:01
2_16			1(12)			2(2,14)				1(6)	4(8 ,9, 11, 12)
HLA-DQA1*05:01
13–27			3(1,5,6)	1(10)	1(12)	3(3,10,12)	1(8)			1(7)	3(1,5,11)
HLA-DQB1*0602
2_16			1(12)			2(2,14)				1(6)	4(8 ,9, 11, 12)
HLA-DQB1*0602
17–31							1(4)
HLA-DQB1*0602
102–116			2(1,2)	2(6,15)	1(8)	4(6, 8, 14, 15)				1(3)	2(1,7)
HLA-DQB1*0301		3(1, 2, 3)	2(6, 9)	2(1, 2)		3(1, 2, 3)					2(6, 9)
HLA-G: HLA-G*01:01, HLA-G*01:02
4_16			1(10)			1(12)				1(4)	4(6 7 9 10)
8_21			3(6,10,11)			1(8)	1(13)			1(12)	5(2, 3, 5 ,6, 10)
22–33				2(1,10)	1(3)	4(1,3,9,10)					1(2)
66–79		1(11)	3(4,5,9)	2(3,12)	1(10)	5(3, 7, 10, 11, 12)	1(6)	1(7)		1(6)	3(2,4,9)
100–112		1(8)	2(2,5)			2(8,11)	1(11)			1(6)	1(2)
HLA-DRB1*0301
96–104		3(1,2,3)	2(6,9)		2(1,2)	3(1,2,3)					2(6,9)
HLA-DRB1*0305
96–104		3(1,2,3)	2(6,9)		2(1,2)	3(1,2,3)					2(6,9)
HLA-DRB1*0309
96–104		3(1,2,3)	2(6,9)		2(1,2)	3(1,2,3)					2(6,9)
HLA-DRB1*11
34–42		2(3,4)	2(7,8)	1(2)	3(2,3,5)	4(2,3,4,5)			1(7)		2(7,8)

The position of each epitope, the number of cleave sites in each epitope as well as the position of the cleave sites in each epitope are mentioned in the table

The results of predicted cleavage positions for 12 commopan class="Chemical">n proteases: B-cell, T-cell, and CTL predicted epitopes

The position of each epitope, the pan class="Chemical">number of cleave sites in each epitope as well as the position of the cleave sites in each epitope are mentioned in the table

Discussion

Although emerging bacterial n class="Disease">viral diseases have caused great catastrophes inclass="Chemical">pan> n class="Species">human history which can affect from a small and localized group to millions of n class="Species">people across continents, several vaccine and therapies have been introduced to control them (Dehghani et al. 2014, 2013).

Fishman et al. (2009) using multivariable logistic regression models ,found 10 n class="Species">HCV-core genpan>e polymorphisms extenpan>sively associated with increased HCC risk (36G/C, 209A, 271C/U, 309A/C, 384U, 408U, 435A/C, 465U, 481A, 546A/C) anclass="Chemical">pan>d onpan>e signpan>ificanpan>tly linked with decreased HCC risk (78U). Menpan>tionpan>ed mutationpan>s related to chanpan>ge in domain1 amino acid sepan> class="Chemical">quence: N11S/T, K12Silent/N, A25V, G69S, Q70R, M91L, and L102P. All current amino acid changes decreased HCC risk except A25V (78U) (Fishman et al. 2009).

In selected sen class="Chemical">quenclass="Chemical">pan>ces, we found that amino acid 11 was T in all sepan> class="Chemical">quences except one sequence in 2013 (P) and in 1a (N).

Positions 12(K) did not show any change; amino acid 25 was P and amino acid 69 was R in all sen class="Chemical">quenpan>ces. Amino acid 70 in nearly all sepan> class="Chemical">quences was R except in 2 sequences in 2006 (Q), 5 sequences in 2013 (3(Q), 2(H), and 1(P)), 4 sequences in 2014 (Q), in 2016 [25(Q), 1(H) and 2 reference sequences (Q)].

Amino acid 91 in most sen class="Chemical">quenpan>ces was C except 2 seclass="Chemical">pan> class="Chemical">quences in 2006(M), 17 sequences in 2013 (13(M), 3(L), 1(F)),50 sequences in 2016(43(M), 7(L)) and three ref sequences, (1b(M), 2a and 5a (L)). In amino acid 102 we did not find any change but in one of the reference sequences, we found one change (5a (G to S)).

Akuta et al. (2007) employed PCR for detecting substitutions of aa 70 and aa 91 in n class="Species">HCV-core genpan>e of genpan>otype 1b by using the mutationpan>-specific primer as anclass="Chemical">pan> important predictor of pan> class="Disease">hepatocarcinogenesis. For wild samples, aa 70 was arginine (R) and aa 91 was n class="Chemical">leucine (L) but for mutant aa 70 was glutamine(Q)/histidine (H) and aa 91 was methionine (M) (Akuta et al. 2007).

Also, Furui et al. (2011) identified aa 70 apan class="Chemical">nd aa 91 substitutions among Japanese volunteer blood donors (Furui et al. 2011).

In terms of aa 70 substitutions, we recognized 2 n class="Chemical">glutamine substitutionpan>s in 2006 seclass="Chemical">pan> class="Chemical">quences and 3 glutamine, 2 histidines, one proline in 2013 sequences; also 4 glutamine substitutions in 2014 sequences, and 25 glutamine, and 1 histidine in 2016. We found several methionine residues in aa 91 in 2006, 2014, and 2016.

Ogata et al. (2002) compared sen class="Chemical">quences of the core protein of Subtype 1b pan> class="Species">HCV strains obtained from patients with and without HCC and found some amino acid mutation sites (Ogata et al. 2002). K23Q, n class="Chemical">Q70R, and T110M substitutions were found by Ogata et al. (2002). In comparison with our results, in all sequences, aa 23 was K, in the majority of sequences aa 70 was R and in 34 sequences it was Q. In aa 110 we did not find any methionine (Ogata et al. 2002).

Akuta et al. (2010) confirmed the role of n class="Chemical">Gln70 (or class="Chemical">pan> class="Chemical">His70) in the efficacy of triple therapy and sustained a virological response, the patient with both genotype non-TT and Gln70 (His70) had the worst sustained virological response. Also, Akuta et al. (2011) by following up twenty-six patients determined the role of Gln70 (His70) substitution in the development of HCC They suggested detection of aa substitutions in the core region before antiviral therapy (Akuta et al. 2010, 2011).

Alestig et al. (2011) showed substitution ipan class="Chemical">n aa 70 of the core was related to treatment response, but that was less important than n class="Gene">IL28B polymorphism. Inpan> our study, 37 sepan> class="Chemical">quences had Gln70 (or His70) substitution (Alestig et al. 2011).

Tokita et al. (2000) approved the role of n class="Species">HCV-core region (pan> class="Chemical">Thr49Pro) to reduce the fluorescence enzyme immunoassay (FEIA) sensitivity. We found 2 sequences in 2013 and one in 2016 with T to P substitutions (Tokita et al. 2000).

Findings of Horie et al. (1999) indicated that alteration from n class="Chemical">glycine to pan> class="Chemical">serine at core codon 45 was dominant in noncancerous liver portions rather than in cancerous liver portions and sera from HCC patients (Horie et al. 1999). Our results did not show any glycine to serine mutation in aa 45 and in all sequences, it was glycine.

Idrees and Ashfan class="Chemical">q (2013) by usinclass="Chemical">pan>g molecular docking software, reported interactionpan>s of amino acid residues pan> class="Chemical">arginine 149, arginine 39, arginine 74, and arginine 78 in HCV-core protein and n class="Chemical">Leu44, Ala71, Ser76 and Pro97 in CXCL6. This finding clues to understanding HCV pathogenesis. Any change in these positions can relate to HCV Infection and HCC. Our results showed no alteration in aa 39 and 74, and in aa 78 just one sequence was glutamine to arginine (Idrees and Ashfaq 2013).

Using a combination of predicted B-cell antibody epitopes by all methods on the immune epitope website and also considering bcepred and ABCpred prediction ,we could define n class="Chemical">three major epitopes (4–20, 50–60, pan> class="Gene">and 100–112) for domain 1.

Ferroni et al. (1993) by using the algorithm of Jameson and Wolf identified four epitopes in n class="Species">HCV-core protein (7–21, 31–45, 49–63, and 99–113).

Harase et al. (1995) analyzed the response to n class="Species">HCV-core protein inclass="Chemical">pan> n class="Species">mice and found a major B cell epitope (21–40).

Pirisi et al. (1995) analyzed n class="Chemical">Sera from 97 pan> class="Species">HCV infected patients and found three 15-mer peptides as antigens in an enzyme immunoassay. They concluded that anti-R15P (50–64, RKTSERSQPRGRRQP) as a potent antigen might help to identify a subgroup at higher risk to develop HCC (Pirisi et al. 1995).

Also, a study on n class="Species">HCV positive blood donclass="Chemical">pan>or by Lechmann et al. (1996) determined a regionpan> (aa 1–24) of domain1 that bounpan>d the anpan>tibodies from the sera of all patients and showed a great potential for detection of HCV infection by using serological B-Cell responses tests.

Comparison of our results with previous studies revealed that all predicted epitopes ipan class="Chemical">n our research have a good potential for future studies of the immune response against n class="Disease">HCV infection, anpan>d are useful for recognitionpan> of all kinds of pan> class="Disease">HCV infections.

Gededzha et al. (2014) used several bioinformatics tools to predict n class="Gene">HLA class I anclass="Chemical">pan>d n class="Gene">HLA class II in n class="Species">HCV genotype 5a. They found three T-cell epitopes of NS3, NS4B, and NS5B.

Some n class="Gene">HLA class II alleles were found in Iranpan>ianpan> pan> class="Species">patients by Samimi-Rad et al. (2015); DRB1*0301, DQA1*0501, DQB1*0201, DRB1*1101, and n class="Gene">DQB1*0301 were demonstrated in patients with HCV clearance, and DRB1*0701, DQA1*0201, DQB1*0602, DRB1*0301, DRB1*11, and DQB1*0201 occurred more frequently in chronic patients (Samimi-Rad et al. 2015).

Khorrami et al. (2015) found a relationship between n class="Gene">HLA-G, pan> class="Gene">IL-10, and response to combined therapy in HCV positive patients. They concluded that HLA-G, and IL-10 have a significant role in response to therapy with IFN-α2α and ribavirin.

Also, n class="CellLine">HLA-A01pan> and HLA-B38 were determined as important alleles associated with Peg-IFN plus n class="Chemical">ribavirin therapy in Egyptian patients by Farag et al. (2013).

Pourhassapan class="Chemical">n et al. (2014) found several HLA alleles associated with HCV in Iranian patients (2011–2013). A2, A3, B35, B38, BW4, CW4and CW7 were the most frequent alleles found by this group.

In accordance with the above-mentioned studies, we collected n class="Gene">HLA alleles associated with pan> class="Disease">HCV infection, and by employing in-silico analysis we established numerous T-cell epitopes for domain1 that can be helpful for future studies to design effective vaccine against HCV genotype 1a, and can provide benefit data for better understanding the role of domain1 in immune response.

Many researchers have proved broader CTL responses to n class="Disease">HCV infection anpan>d the usefulnclass="Chemical">pan>ess of CTL epitopes mapping to develop therapeutic intervenpan>tionpan>s or vaccines (Sabet et al. 2014; Saeedi et al. 2014; Jazayeri anpan>d Carmanpan> 2005; Arashkia et al. 2011). Inpan> our research, we utilized reliable software to predict CTL epitopes anpan>d extracted data that canpan> be useful for vaccine developmenpan>t studies.

By considering results of phosphorylation sites prediction by 3 programs, we concluded that 3 sites (15, 49, n class="Gene">and 116) were the main phosphorylationpan> sites in domain 1. Aminclass="Chemical">pan>o acid 116 is a n class="Chemical">serine that located in n class="Chemical">Arg-Arg-Arg-Ser-Arg region; this region was similar to the usual target sequence for protein kinase A [Arg-Arg-X-Ser/Thr-X]. Two threonine amino acid residues (15, 49) were calculated as protein kinase c target sites where this kinase acts through the phosphorylation of hydroxyl groups of amino acid residue. Previous studies indicated that core protein is phosphorylated by PKA and PKC.

Core phosphorylation regulates the suppressive activity of n class="Species">HCV-core proteinclass="Chemical">pan> onpan> HBV genpan>e replicationpan> anpan>d expressionpan>. Also, it has beenpan> shownpan> that phosphorylationpan> in core relates to the nuclear localizationpan> of the core protein. They demonpan>strated three serine residues (Ser-53, Ser-116, and Ser-99) as the potential phosphorylated sites in core protein, that were similar to NetPhos software results in our research (Yassin 2001; Shih et al. 1995; Lu and Ou 2002).

Secondary structure predictiopan class="Chemical">n indicated that the majority of domain1 was the random coil, and all B-cell epitopes and important mutations placed on random coil structure.

Tertiary structures were designed by n class="Chemical">three signclass="Chemical">pan>ificant anpan>d reliable onpan>line programs, but just onpan>e of them provided a reliable anpan>d high-pan> class="Chemical">quality protein structure model for domain1. The quality and reliability of models were confirmed by QMEAN and RAMPAGE software.

By examining all the predictions core epitopes with and without signal peptide, we found out that there is no difference between these two different strategies of analysis (Pene et al. 2009; Tn class="Chemical">argett-Adams et al. 2008; Ma et al. 2007; Okamoto et al. 2008; Oehler et al. 2012).

Based on previous studies core proteipan class="Chemical">n has a C terminal signal peptide (170–191) and because the focus of our study was on domain 1, it was expected that the deletion of this region could not affect the epitope perdition results.

Digestion apan class="Chemical">nalysis to predict possible proteases was shown that each epitope can be digested by at least 5 selected proteases which can have a significant effect on the reduction of the half-life of epitopes.

Conclusion

Finally our investigations in this research provided comprehensive data about fren class="Chemical">quenpan>t mutationpan>s in domain1, anclass="Chemical">pan>d as a first report can be useful for future study about signpan>ificanpan>t mutationpan>s anpan>d their role in therapeutic pathway anpan>d responpan>se to anpan>tiviral therapy for Iranpan>ianpan> pan> class="Species">patients. Also, identification of domain1 properties provides practical information for domain1 cloning and more researches.