A bioinformatic approach to investigating cytokine genes and their receptor variants in relation to COVID-19 progression.

Severe acute respiratory syndrome coronavirus 2 infection produces a wide spectrum of manifestations, ranging from no symptom to viral pneumonia. This study aimed to determine the genetic variations in cytokines and their receptors in relation to COVID-19 pathogenesis using bioinformatic tools. Single nucleotide polymorphisms (SNPs) of genes encoding the cytokines and cytokine receptors elevated in patients with COVID-19 were determined from the National Biotechnology Information Center website (using the dbSNP database). Missense variants were found in 3 cytokine genes and 10 cytokine receptor genes. Computational analyses were conducted to detect the effects of these missense SNPs via cloud-based software tools. Also, the miRSNP database was used to explore whether SNPs in the 3'-UTR altered the miRNA binding efficiency for genes of cytokines and their receptors. Our in silico studies revealed that one SNP in the vascular endothelial growth factor receptor 2 (VEGFR2) gene was predicted as deleterious using sorting intolerant from tolerant. Also, the stability of VEGFR2 decreased in the I-Mutant2.0 (biotool for predicting stability changes upon mutation from the protein sequence or structure) prediction. It was suggested that the decrease in VEGFR2 function (due to the rs1870377 polymorphism) may be correlated with the progression of COVID-19 or contribute to the pathogenesis. Moreover, 27 SNPs were determined to affect miRNA binding for the genes of cytokine receptors. CXCR2 rs1126579, TNFRSF1B rs1061624 and IL10RB rs8178562 SNPs were predicted to break the miRNA-mRNA binding sites for miR-516a-3, miR-720 and miR-328, respectively. These miRNAs play an important role in immune regulation and lung damage repair. Further studies are needed to evaluate the importance of these miRNAs and the SNPs.

INTRODUCTION

Severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), which belongs to the Coronaviridae family, causes respiratory and gastrointestinal infections. The World Health Organization named the disease caused by this virus COVID‐19, which is an acronym for ‘coronavirus disease 2019’, while the agent was named SARS‐CoV‐2 due to its similarity to SARS‐CoV (Bassetti et al., 2020). Since the science and medical community has not faced such a widespread epidemic before, local experiences come to the forefront of managing this situation (Rombolà et al., 2020). Patients with severe symptoms who require hospitalization for SARS‐CoV‐2 infection include men, old people, smokers, patients with obesity and those with common comorbidities such as cardiovascular diseases, diabetes and chronic lung disease (Yang et al., 2020). However, to reduce the mortality rate of COVID‐19, further investigation is still needed to find effective indicators for assessing the severity and clinical progression of the disease. Some of the patients show only mild fever, cough or muscle soreness, while some patients' conditions deteriorate in the later stages and result in death due to acute respiratory distress syndrome (ARDS) and multiple organ failure (Guo et al., 2020). Huang et al. (2020) reported the clinical features and cytokine profile of patients with COVID‐19 in Wuhan, China, and suggested that a cytokine storm could be associated with the severity of the disease. In addition, Xu et al. (2020) examined biopsy samples from the deceased, and interstitial mononuclear inflammatory infiltrates predominated by lymphocytes were seen in both lungs. The SARS‐CoV‐2 infection causes a sequential release of specific cytokines that cause significant damage to the pulmonary epithelium, resulting in ARDS, sepsis and organ failure (Mehta et al., 2020).

In Huang et al.'s (2020) study, initial plasma IL1RA, IL1B, IL7, IL8, IL9, IL10, basic FGF, GCSF, GMCSF, IFNγ, IP10, MCP1, MIP1A, MIP1B, PDGF, TNFα and vascular endothelial growth factor (VEGF) concentrations were higher in patients with COVID‐19 than in healthy controls. In addition, of a total of 81,385 cases of COVID‐19 reported by the Chinese Center for Disease Control and Prevention, 81% were mild, 14% severe and 5% critical (Wu & McGoogan, 2020). Thus, genetic variations in cytokines and their receptors could play an important role in the progression or severity of COVID‐19 infection.

The vascular structure in the respiratory system plays an important role in maintaining the physiological functions of expansive capacity and major plasticity. Various pathologic conditions, including COVID‐19, increase the permeability of vascular endothelial cells, expression of adhesion molecules, migration and proliferation of endothelial cells, and infiltration of inflammatory cells (McDonald, 2001). VEGF is considered the most important factor (Riedel et al., 2002) because of the increase in the inflammatory process and serum levels of VEGF in patients with COVID‐19. VEGF produces this effect by binding to VEGF receptor type 1 (VEGFR1) or VEGF receptor type 2 (VEGFR2), which have tyrosine kinase activity (Shibuya & Claesson‐Welsh, 2006). VEGF‐R2 is regarded as the main signalling receptor for VEGF bioactivity (angiogenesis, proliferation and permeability) and can cause proliferation in cells lacking VEGFR1 (Carmeliet et al., 2001). Downstream signal transduction pathways are triggered by VEGFR2 receptor kinase activity—which promotes the proliferation, migration and differentiation of endothelial cells and enhances the permeability of the microvasculature. Alveolar apoptosis and emphysema occur when VEGF activity is inhibited by VEGFR2 in rats (Kasahara et al., 2000).

Single nucleotide polymorphisms (SNPs) can be located in other gene regions—such as 5′‐ or 3′‐UTRs, introns or promoters and the exonic region. Genetic variations in 3′‐UTRs can modify gene expression via miRNA binding, protein–mRNA interactions, gene expression disruption and polyadenylation; therefore, SNPs in 3′‐UTRs are very important and arouse the interest of researchers. Furthermore, it has been shown that SNPs in 3′‐UTRs can affect miRNA functions by changing thermodynamic properties of the hybridization site and the secondary structure of 3′‐UTRs, lowering binding yield, exchanging miRNA recognition elements and, probably, creating new binding sites or enhancing binding efficiency between the target site and miRNA (Schwerk & Savan, 2015; Steri et al., 2018). These 3′‐UTR‐located SNPs have been found to be useful tools for the development of medicine, assessment of disease susceptibility and monitoring of the clinical symptoms of patients in several studies (Ding et al., 2018).

There are a significant number of SNPs in genes encoding cytokines that are high in patients with COVID‐19, and the process of verifying the potential relationship between SNPs and diseases in the laboratory is costly and, most importantly, time‐consuming. In silico analyses allow for narrowing the regions of potential SNP targets for experimental validation. Using bioinformatic tools, the present study aimed to determine the genetic variations in cytokines and cytokine receptors that are possibly related to COVID‐19 pathogenesis.

MATERIALS AND METHODS

The SNPs for genes coding the cytokines and their receptors that were elevated in patients with COVID‐19 were chosen from the dbSNP database, which is available on the National Biotechnology Information Center website (http://www.ncbi.nlm.nih.gov/SNP). Variant analysis was carried out for the SNPs (synonymous and nonsynonymous) in the coding region and the untranslated regions with MAF > 0.15. We analysed the missense SNPs in 3 cytokine genes and 10 cytokine receptor genes using sorting intolerant from tolerant (SIFT) to predict the deleterious and tolerated SNPs (https://sift.bii.a‐star.edu.sg/www/SIFT_dbSNP.html). SIFT uses sequence homology or physical properties to predict the effects of amino acid substitution on protein function and, hence, potential alteration on phenotype (Kumar et al., 2009). Further analysis was conducted for the deleterious SNPs identified with SIFT by PolyPhen and I‐Mutant2.0 database. PolyPhen prediction is based on a series of features—including phylogenetic, structural and sequence annotation information characterizing a substitution. PolyPhen classifies the SNPs as ‘probably damaging’, ‘possibly damaging’ or ‘benign’ (Ramensky et al., 2002). For the prediction of the missense SNP impact on the stability of the protein, the I‐Mutant2.0 database was used (Capriotti et al., 2005). Furthermore, miRNA binding efficiency was also affected by the SNPs in the 3′‐UTR. Therefore, miRSNP was used to predict whether the SNPs affected the miRNA binding efficiency for the cytokine genes and receptors.

RESULTS

The number of the synonymous and nonsynonymous coding SNPs and the SNPs from the untranslated regions of the genes coding cytokines and their receptors that were possibly related to COVID‐19 pathogenesis is listed in Table 1. SIFT was used for the functional significance analysis of the SNPs. The prediction results of the missense SNPs by SIFT are presented in Table 2. One SNP in the VEGFR2 gene was predicted as deleterious, and 13 SNPs were predicted as ‘tolerable’ by SIFT. The other SNP tools such as PolyPhen and I‐Mutant2.0 were used for further analysis of this deleterious SNP. PolyPhen predicted that as benign, where the I‐Mutant2.0 prediction showed decreased stability (Table 2).

TABLE 1

Number of the synonymous and nonsynonymous coding SNPs and the SNPs from the untranslated regions and intron of the genes coding cytokines and cytokine receptors that were high plasma concentrations in COVID‐19 patients

Gene name	3′‐UTR	5′‐UTR	Intron upstream and downstream transcript variant	Synonymous variant	Missense variant
Cytokine genes
CCL2	1	—	2	1	—
CCL3	2	—	7	1	—
CCL4	—	—	29	1	2
CSF2	—	—	8	—	1
CSF3	2	—	7	1	—
CXCL8	2	—	6	—	—
CXCL10	3	—	17	—	—
IFNG	—	—	9	—	—
IL1B	—	1	10	—	—
IL1RN	2	4	143	2	—
IL7	15	—	109	—	—
IL9	—	—	6	—	—
IL10	1	1	12	—	—
PDGFB	—	1	63	—	—
TNF	—	—	2	—	—
VEGFA	2	1	40	—	—
Cytokine receptor genes
CCR1	—	—	4	—	—
CCR2	3	1	15	—	1
CCR4	2	—	5	—	—
CCR5	2	2	14	—	—
CSF2RA	9	1	420	1	—
CSF3R	—	1	22	2	—
CXCR1	—	—	1	—	—
CXCR2	3	1	15	1	—
CXCR3	—	—	2	—	—
FLT1	16	—	254	2	—
IFNGR1	—	3	16	—	—
IFNGR2	2	3	59	—	1
IL1R1	3	7	260	—	—
IL1R2	1	17	110	—	—
IL7R	11	—	53	—	3
IL9R	—	—	30	—	—
IL10RA	1	—	17	—	1
IL10RB	4	4	85	—	1
KDR	—	1	89	—	2
TNFRSF1A	1	1	22	1	—
TNFRSF1B	5	—	47	1	1

TABLE 2

SIFT, PolyPhen‐2 and I‐Mutant‐2.0 results of missense SNPs of genes encoding the cytokines and cytokine receptors that was elevated in COVID‐19 patients

Gene name	Gene ID	SNP	Allele change	Amino acid change	SIFT prediction	Polyphen‐2 prediction	I‐Mutant‐2.0 prediction
Cytokine genes
CCL4	ENSG00000129277	rs1049807	A/G	N41S E79E	Tolerated
CCL4	ENSG00000129277	rs1719152	T/A	S80T N41K	Tolerated
CSF2	ENSG00000164400	rs25882	T/C	I117T	Tolerated
Cytokine receptor genes
CCR2	ENSG00000121807	rs1799864	G/A	V64I	Tolerated
IFNGR2	ENSG00000159128	rs9808753	A/G	Q83R, Q64R	Tolerated
IL7R	ENSG00000168685	rs1494558	T/C	I66T	Tolerated
IL7R	ENSG00000168685	rs6897932	C/T	T244I	Tolerated
IL7R	ENSG00000168685	rs1494555	G/A	V138I	Tolerated
IL10RA	ENSG00000110324	rs2229113	A/G	R351G R331G R202G	Tolerated
TNFRSF1B	ENSG00000028137	rs1061622	T/G	M196R	Tolerated
IL10RB	ENSG00000243646	rs2834167	A/G	K47E	Tolerated
KDR	ENSG00000128052	rs1870377	T/A	Q472H	Deleterious	Benign	Decrease stability
KDR	ENSG00000128052	rs2305948	C/T	V297I	Tolerated

The possible alterations caused by the 3′‐UTR SNPs were investigated for the miRNA binding efficiency in the listed genes that were suggested to have potential roles in COVID‐19. We predicted that the 27 SNPs affected the miRNA binding for the cytokine receptor genes and 10 SNPs for cytokine genes by in silico analysis (Table 3). By using this software, one can predict SNPs' effect on miRNA binding sites, which can be decreased, enhanced, created or broken miRNA binding.

TABLE 3

miRSNP results of SNPs' miRNA binding efficiency for cytokine and cytokine receptor genes

Gene	SNP	Allele	Effect
			Decrease	Enhance	Create	Break
Cytokine genes
CCL2	rs13900	C	hsa‐miR‐3163	hsa‐miR‐374a‐5p hsa‐miR‐374b‐5p	hsa‐miR‐4761‐5p hsa‐miR‐624‐3p
CCL3	rs8951	G	hsa‐miR‐5002‐3p	hsa‐miR‐4672	hsa‐miR‐3929 hsa‐miR‐4419b hsa‐miR‐4438 hsa‐miR‐4478 hsa‐miR‐4502
	rs1063340	G			hsa‐miR‐3179 hsa‐miR‐3202 hsa‐miR‐4716‐3p hsa‐miR‐4723‐5p hsa‐miR‐4747‐5p hsa‐miR‐5196‐5p hsa‐miR‐5698	hsa‐miR‐1292
CSF3	rs2827	C	hsa‐miR‐3653	hsa‐miR‐3658		hsa‐miR‐548ad
	rs1042658	C		hsa‐miR‐2355‐5p	hsa‐miR‐5586‐5p	hsa‐miR‐1247‐5p
IL8	rs1126647	A		hsa‐miR‐944
CXCL10	rs3921	C	hsa‐miR‐5002‐5p			hsa‐miR‐509‐3p hsa‐miR‐591
	rs34836828	Deletion		hsa‐miR‐145‐3p
VEGFA	rs3025040	C			hsa‐miR‐199a‐5p hsa‐miR‐199b‐5p hsa‐miR‐4676‐5p hsa‐miR‐575
	rs10434	A		hsa‐miR‐4727‐5p	hsa‐miR‐3677‐5p	hsa‐miR‐3545‐5p hsa‐miR‐5693 hsa‐miR‐660‐3p
Cytokine receptor genes
CCR2	rs743660	A		hsa‐miR‐4786‐3p
CCR5	rs746492	T		hsa‐miR‐5007‐3p	hsa‐miR‐4524a‐3p hsa‐miR‐589‐3p	hsa‐miR‐3133
CXCR2	rs1126579	C		hsa‐miR‐5193	hsa‐miR‐138–1‐3p	hsa‐miR‐516a‐3p hsa‐miR‐516b‐3p
	rs1126580	A	hsa‐miR‐4524b‐3p		hsa‐miR‐5096
FLT1	rs2296283	C			hsa‐miR‐3135b hsa‐miR‐3940‐3p	hsa‐miR‐1538 hsa‐miR‐4731‐5p hsa‐miR‐4745‐3p
	rs2296284	C	hsa‐miR‐3943 hsa‐miR‐4313		hsa‐miR‐4293	hsa‐miR‐1234
	rs3209052	A				hsa‐miR‐4789‐3p hsa‐miR‐582‐5p
	rs3751397	T	hsa‐miR‐3662			hsa‐miR‐548a‐3p hsa‐miR‐548ar‐3p hsa‐miR‐548e hsa‐miR‐548f
	rs7326277	G			hsa‐miR‐193b‐5p hsa‐miR‐4446‐5p	hsa‐miR‐193a‐5p
	rs7337610	G		hsa‐miR‐589‐3p	hsa‐miR‐448
	rs9551465	T	hsa‐miR‐223‐5p
	rs17086617	G				hsa‐miR‐224‐5p
	rs35779457	Deletion	hsa‐miR‐876‐3p			hsa‐miR‐4495
	rs55875014	G		hsa‐miR‐3916 hsa‐miR‐5197‐3p		hsa‐miR‐3065‐5p hsa‐miR‐3529‐3p hsa‐miR‐3928
	rs56340749	Deletion	hsa‐miR‐4446‐3p hsa‐miR‐4498	hsa‐miR‐194‐3p hsa‐miR‐5001‐5p		hsa‐miR‐1225‐5p
	rs56791288	Deletion		hsa‐miR‐4531
IFNGR2	rs1059293	C		hsa‐miR‐493‐5p
IL1R1	rs2110726	G			hsa‐miR‐4534 hsa‐miR‐4802‐5p
	rs3732131	G			hsa‐miR‐4762‐3p
	rs3917324	C		hsa‐miR‐4716‐5p hsa‐miR‐4776‐3p hsa‐miR‐4781‐3p	hsa‐miR‐191‐3p hsa‐miR‐604 hsa‐miR‐647	hsa‐miR‐1587 hsa‐miR‐378g hsa‐miR‐4492 hsa‐miR‐4498 hsa‐miR‐4505 hsa‐miR‐5001‐5p hsa‐miR‐762
IL10RB	rs1058867	A	hsa‐miR‐219–1‐3p			hsa‐miR‐377‐5p
	rs3171425	A		hsa‐miR‐328	hsa‐miR‐1282 hsa‐miR‐4655‐3p
	rs8178562	A		hsa‐miR‐4252 hsa‐miR‐5008‐3p		hsa‐miR‐218‐5p hsa‐miR‐328 hsa‐miR‐636
TNFRSF1B	rs3397	C	hsa‐miR‐3126‐5p	hsa‐miR‐5581‐5p	hsa‐miR‐329 hsa‐miR‐362‐3p	hsa‐miR‐122‐3p
	rs1061624	A/C/G		hsa‐miR‐3188	hsa‐miR‐5003‐5p
		C, T→A T→G		hsa‐miR‐3692‐3p
		G→C,A	hsa‐miR‐3692‐3p
		C	hsa‐miR‐922			hsa‐miR‐523‐3p
		A		hsa‐miR‐922	hsa‐miR‐523‐3p	hsa‐miR‐639 hsa‐miR‐720
	rs1061628	A			hsa‐miR‐4715‐5p	hsa‐miR‐3680‐3p
		C				hsa‐miR‐4715‐5p
	rs5746065	A	hsa‐miR‐4786‐3p	hsa‐miR‐1299 hsa‐miR‐671‐5p	hsa‐miR‐4731‐5p hsa‐miR‐486‐3p	hsa‐miR‐5589‐5p

DISCUSSION

We determined the genetic variations of genes coding cytokines and receptors in relation to COVID‐19 by using bioinformatic tools. There are four missense SNPs in genes encoding cytokines that had high plasma concentrations in patients with COVID‐19. But SIFT analysis predicted that these variations are tolerable and not expected to affect the protein function. However, of the 10 missense polymorphisms found in genes coding the receptors, VEGFR2 gene Q472H (rs1870377) polymorphism was predicted to have a deleterious effect by SIFT and I‐Mutant2.0 prediction database predicted that this polymorphism decreased the stability of the protein.

The VEGFR2 gene consists of 26 exons, is located in 4q11–q12 and encodes 1,356 amino acids. Missense substitution (c.1416A > T) causes Q472H change in the fifth extracellular Ig‐like motifs (Glubb et al., 2011). VEGF has an important function in suppressing the apoptosis cascade and reducing oedema formation by decreasing the increased endothelial permeability following the intratracheal application of inflammatory stimuli. Koh et al. (2007) reported that VEGF is a major protective factor for the damaged lung during the progression of ARDS. The decrease in VEGFR2 function due to the rs1870377 polymorphism may be the reason for vascular dysfunction—including impaired endothelial cell survival, endothelial cell damage and abnormal vascular repair, contributing to the progression of COVID‐19 and the pathogenesis. However, there are no sequence data for rs1870377 variants from patients with COVID‐19 and the frequency of this variant is close to 50% in East Asian, Vietnamese and Korean populations. This limits the impact of our findings.

Gene regulation has an essential role in host defence against pathogens, and its dysregulation has been demonstrated in different infectious diseases or disease progression (Chandan et al., 2020). 3′‐UTR polymorphism in genes encoding cytokines or their receptors was found in higher levels in patients with COVID‐19 (Table 3). The effect of SNPs on miRNA binding efficiency and the variant that is responsible for the described effect is shown in Table 3. Among these polymorphisms, CXCR2 rs1126579, TNFRSF1B rs1061624 and IL10RB rs8178562 are particularly remarkable because these SNPs would break the miRNA‐mRNA binding sites for miR‐516a‐3p, miR‐720 and miR‐328, respectively. It has been suggested that the main cause of lung injury during a response to SARS‐CoV‐2 is an increase in these pro‐inflammatory cytokines and the dysregulation of the immune response.

Narożna et al. (2017) reported that miR‐328 represents a potent modifier of the complex process of wound repair in bronchial epithelial cells and inhibition of miR‐328 interrupts the repair process. In addition, Wu et al. (2019) demonstrated that miR‐516a‐3p expression knockdown could inhibit cell proliferation, invasion, migration and wound repair but promote apoptosis in lung adenocarcinoma cells. Also, previous studies have reported that miR‐328 plays important role in regulating the expression of genes associated with cell–cell interactions, transport across the membranes (Li et al., 2011), migration and cell adhesion (Ishimoto et al., 2014), and calcium‐dependent processes such as cell division, cell motility and cell death (Lu et al., 2010). Hence, the increased expression of these miRNAs is inevitable to repair lung damage that is associated with COVID‐19. We can speculate that high expression of miR‐516a‐3p and miR‐328, as a result of the repair of the lung damage, by reducing the cytokine levels, might be obstructed by the related SNPs. However, further studies are needed to determine the role of these miRNAs in COVID‐19 pathogenesis.

The total number of NK and CD8+ T cells decreased markedly in patients with SARS‐CoV‐2 infection (Zheng et al., 2020). miR‐720 regulates TCR‐mediated proliferation of primary human CD8+ T cells and has an important role in immune regulation (Wang et al., 2015). The upregulation of miR‐720 in CD8+ T cells may play a role in the progression of COVID‐19. It was also found in our study that rs1061624 polymorphism in the TNFRSF1B gene disturbed the miR‐720 regulator effect on TNFRSF1B expression. Therefore, this is a valid candidate for further study on the pathogenesis and progression of COVID‐19.

In conclusion, our analysis, using the bioinformatic approach, showed that VEGFR2 rs1870377 polymorphism comes into prominence according to SIFT and the I‐Mutant2.0 database. Also, CXCR2 rs1126579, TNFRSF1B rs1061624 and IL10RB rs8178562 attracted attention because it was predicted that these SNPs could break the miRNA‐mRNA binding sites for miR‐516a‐3, miR‐720 and miR‐328, which are important miRNAs in immune regulation and repair of damage in the lungs.

CONFLICT OF INTEREST

The authors report no conflicts of interest in this work.

AUTHOR CONTRIBUTIONS

Güneş Çakmak Genç carried out the literature review, while Sevim Karakaş Çelik conducted the bioinformatic analyses. Güneş Çakmak Genç, Sevim Karakaş Çelik and Ahmet Dursun performed the evaluation and discussion of the results.