Bioinformatics analysis reveals molecular connections between non-alcoholic fatty liver disease (NAFLD) and COVID-19.

The ongoing coronavirus disease 2019 (COVID-19) pandemic caused by SARS-CoV-2 has devastatingly impacted people's lives. Non-alcoholic fatty liver disease (NAFLD) is fatal comorbidity of COVID-19 seen with potential risk factors to develop severe symptoms. This research focuses on determining and elucidating the molecular factors and connections that might contribute to the severity of SARS-CoV-2 infection in NAFLD patients. Here, we comprehensively inspected the genes involved in NAFLD and SARS-CoV-2 entry factors (SCEFs) found by searching through the DisGeNet database and literature review, respectively. Further, we identified the SCEFs-related proteins through protein-protein interaction (PPI) network construction, MCODE, and Cytohubba. Next, the shared genes involved in NAFLD and SARS-CoV-2 entry, and hub gene were determined, followed by the GO and KEGG pathways analysis. X2K database was used to construct the upstream regulatory network of hub genes, as well as to identify the top ten candidates of transcription factors (TFs) and protein kinases (PKs). PPI analysis identified connections between 4 top SCEFs, including ACE, ADAM17, DPP4, and TMPRSS2 and NAFLD-related genes such as ACE, DPP4, IL-10, TNF, and AKT1. GO and KEGG analysis revealed the top ten biological processes and pathways, including cytokine-mediated signaling, PI3K-Akt, AMPK, and mTOR signaling pathways. The upstream regulatory network revealed that AKT1 and MAPK14 as important PKs and HIF1A and SP1 as important TFs associated with AKT1, IL-10, and TNF. The molecular connections identified between COVID-19 and NAFLD may shed light on discovering the causes of the severity of SARS-CoV-2 infected NAFLD patients.

Introduction

Since December 2019, following several reported cases of pneumonia of unknown cause, a novel coronavirus (2019-nCoV) was confirmed as the causative pathogen on January 7, 2020. The novel coronavirus nominated as Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) emerged in the city of Wuhan, China, and caused a prevalence of unfamiliar viral pneumonia. This novel coronavirus disease, also known as coronavirus disease 2019 (COVID-19), has spread to every corner of the globe and has been blamed for million deaths (Hu et al. 2020). Despite the efforts of scientists and the health care society, many countries are still facing the consequences of the infamous virus (Wu and McGoogan 2020). To this end, different variants of the virus have been detected, and more recently, the Omicron variant was identified in many countries, and it seems that it likely will spread more quickly than the original SARS-CoV-2 virus and other variants like Delta. Among COVID-19 patients, the disease severity was substantially associated with weakened immune response or comorbidities (Fang et al. 2020). Risk factors such as age, sex, lifestyle, and pre-existing conditions, for instance, chronic liver disease, play a significant role in COVID-19 severity (Gao et al. 2021). A single-center case series involving 138 patients showed that comorbidities in COVID-19 patients and old age significantly impacted them, thus suffering a severe form of infection, and they were more likely to be admitted to the intensive care unit (ICU) (Wang et al. 2020).

Several key cell entry mechanisms and factors of SARS-CoV-2 have recently been identified that possibly contribute to the cell infectivity and immune evasion of the virus (Shang et al. 2020; Jackson et al. 2022). However, the underlying mechanisms of molecular action and entry of SARS-CoV-2 are not completely clear yet. It has been shown that SARS-CoV-2 enters into its host cells by interacting with receptors/coreceptors, as well as with other cofactors, through its spike (S) protein that further facilitates fusion between viral and cellular membranes (Peng et al. 2021). Among the suggested factors and mechanisms, Angiotensin-Converting enzyme2 (ACE2) receptor, type II transmembrane serine protease (TMPRSS2), and a disintegrin and metalloproteinase 17 (ADAM17) seem to play a more important role in SARS-CoV-2 infection (Zipeto et al. 2020; Seyedpour et al. 2021; Jackson et al. 2022). ACE2 and TMPRSS2 expressed in different organs, tissues, and specifically in the lungs as a primary receptor of SARS-CoV-2, work in the first line of viral activity and the viral entrance (Baughn et al. 2020). Viral entry in SARS-CoV2- into the host cell occurs while they access densely glycosylated S proteins of ACE2 (Wrapp et al. 2020a, b). By binding to S proteins and cleaving the C-terminal portion of ACE2 by proteases such as TMPRSS2 and FURIN, which are expressed in lung tissue, they access host cells (Pinto et al. 2020). Recent studies showed that disruption of ACE2 by SARS-CoV-2 may manifest as short- or long-term pathophysiological problems in the lungs, brain, heart, gastrointestinal tract, pancreas, and oral cavity (Ashraf et al. 2021).

NAFLD (Non-Alcoholic Fatty Liver Disease) is a multi-organ and multifactorial disease caused by a build-up of fat in the liver and it is usually seen in people who are overweight. It has been reported that NAFLD causes the expression of SARS-CoV2 entry factors in different tissues (Meijnikman et al. 2021). A recent analysis suggested that NAFLD may result in a malfunctioned immune response, making these patients significantly more susceptible to experience a severe stage of COVID-19 than non-NAFLD patients (Ji et al. 2020). Understanding the molecular connection between NAFLD-related genes and the COVID-19 mechanism can help us to facilitate treatment and give us a better clinical explanation about how pre-existing diseases lead to an even worse infection.

Here, we aim to represent the possible reason for the severity of SARS-CoV2 infected patients with pre-existing NAFLD at the molecular level by analyzing the co-expression networks of the SARS-CoV2 entry factors (SCEFs). First, we selected 12 important SCEFs, then identified most important genes/proteins that potentially interact with the selected SCEFs. Next, we analyzed the NAFLD-related and SCEFs-related genes, and their shared genes were determined as NS genes. Finally, the hub genes were identified and analyzed using different bioinformatics tools.

Methods and materials

COVID-19 and NAFLD datasets

To obtain a list of genes involved in SARS-CoV2 entry factors (SCEFs), a comprehensive literature search was performed. COVID-19-related words such as “COVID-19”, “SARS-CoV2 entry factors”, “SARS-CoV2 receptors”, “Disease severity”, “Risk factors” were used as the keywords. In order to find the NAFLD-related genes, the DisGeNet database was used (http://www.disgenet.org/). The DisGeNET database integrates information of human gene-disease associations (GDAs) and variant-disease associations (VDAs) (Piñero et al. 2021). Here, three UMLS IDs, including “C0015695”, “C0400966”, and “C3241937” were searched in the DisGeNET database.

Protein-protein interaction network analysis and selection of NS and hub genes

To find the important genes/proteins interacting with SCEFs, the protein-protein interactions (PPI) network for each of 12 selected SCEFs was analyzed using STRING online database (https://www.string-db.org), and BioGrid online tool (https://thebiogrid.org). Then, Cytoscape V3.8.2 (http://www.cytoscape.org/) was used for the demonstration of the PPI network (An interaction with a combined score > 0.4 was set as the cut-off value). To find the PPI sub-network, the MCODE clustering algorithm was used (The criteria were set as follows: degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max depth = 100). Additionally, Cytohubba was used to single out the top 10 genes/ proteins with more interactions. The factors of the first cluster (found by MCODE) and the top 10 genes/proteins (identified by Cytohubba) were determined as SCEFs-related genes. The common genes between NAFLD-related and SCEFs-related genes were determined as NS genes. Additionally, SCEFs interacting with ≥ 9 NS genes were identified. The selected SCEFs and their related genes were designated as Hub genes.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis

GO and KEGG pathway for the hub genes were analyzed using the EnrichR online tool (https://maayanlab.cloud/Enrichr/). GO enrichment analysis is categorized into three groups: Biological Process (BP), Molecular Function (MF), and Cellular Component (CP). Also, the KEGG pathway was used to classify genes and biological pathways. Adj P-value < 0.05 was set as the cut-off criterion.

Upstream regulatory network analysis

The X2K database (https://maayanlab.cloud/X2K/) was used to identify the transcription factors (TFs), protein kinases (PKs), and their related intermediate proteins that potentially regulate hub genes expression. Cytoscape was used to visualize the gene regulatory network (GRN) for hub genes. The combined Adj P-value of < 0.05 was considered as a significance threshold.

Micro RNA regulatory analysis

MicroRNAs are small non-coding RNAs that control the expression of target genes at the post-transcriptional level. miRTarBase_2017 database from the EnrichR was used to construct the miRNA-target interactions (MTIs) network that potentially regulate the hub genes expression. We selected the top ten miRNAs with Adj P-value < 0.05 and manifested the MTIs network using Cytoscape.

Results

Identification of COVID19 and NAFLD-related genes

The genes functioning as SCEFs and viral replication were investigated at recent studies. Angiotensin-converting enzyme 2 (ACE2), A disintegrin and metalloprotease 17 (ADAM17), Cathepsin L (CTSL), basigin (BSG), C-Type Lectin Domain Family 4 Member G (CLEC4G), C-Type Lectin Domain Family 4 Member M (CLEC4M), Cathepsin B (CTSB), dipeptidyl peptidase 4 (DPP4), Glucose regulatory protein 78 (GRP78) also called Heat Shock Protein Family A (Hsp70) Member 5 (HSPA5), Lymphocyte Antigen 6 Family Member E (LY6E), Zinc Finger CCHC-Type And RNA Binding Motif Containing 1 (ZCRB1) also called MADP1, and Transmembrane Serine Protease 2 (TMPRSS2) were determined as the Top SARS-CoV-2 entry related factors.

Additionally, NAFLD-related genes were obtained using the DisGeNet database. Analysis of “C0015695”, “C0400966”, and “C3241937” revealed that 1561 genes participate in NAFLD development and pathogenesis.

Identifying NS and hub genes

In order to find the important factors that interact with each selected SCEFs, we used STRING and BioGrid online databases. The constructed PPI network, MCODE and Cytohubba revealed that each SCEFs interacted with several specific proteins (Table 1). Among them, ACE2 and TMPRSS2, as two known SCEFs, interacted with 19 and 16 different key proteins, respectively. Additionally, CLEC4M interacted with only 11 proteins, whereas MADP1 had the highest number of interactions with other proteins (31 interactions).

Table 1

SARS-CoV2 entry factors (SCEFs) found by literature review and their important interacting genes/proteins

SCEFs	Interacting genes/proteins
ACE2	AGTR1, AGT, ACE2, KNG1, REN, GCG, EDN1, MME, ACE, AGTR2, XPNPEP2, PRCP, APLN, ENPEP, CMA1, MEP1B, MEP1A, ATP6AP2, DPP4
ADAM17	JAG1, EGF, TNFRSF1A, ERBB2, EGFR, NOTCH1, IL10, JAG2, ADAM17, TNF, IKBKG, IKBKB, MAP3K7, TAB2, TRADD, RIPK1, TRAF2, BIRC2
CTSL	HLA-DQA1, HLA-DPB1, HLA-DRB5, CTSL, HLA-DRA, HLA-DRB1, HLA-DQB1, HLA-DQA2, HLA-DPA1, HLA-DQB2, CTSS, HLA-DMB, HLA-DMA, PTPN22, CD4, CD74
BSG	ITGA5, EGFR, SRC, ITGA3, BSG, ITGB1, ITGA6, ITGA9, TLN1, ITGAV, ITGA4, ITGA8
CLEC4G	LCP2, CLEC5A, CLEC4G, SYK, TREM1, CD300LB, TYROBP, HCST, TREM2, SIRPB1, KLRD1, KLRC2
CLEC4M	RHOA, ARHGAP35, ARHGEF11, CLEC4M, ARHGEF12, ROCK1, CD209, RHOC, ARHGDIA, RTKN, ROCK2
CTSB	HLA-DRB5, HLA-DPB1, HLA-DQA1, HLA-DQB2, HLA-DPA1, HLA-DQB1, HLA-DRB1, CTSS, CTSB, CD74, CTSL, CTSD, HLA-DRA, HLA-DMA, HLA-DMB, CD4
DPP4	INS, AGT, DPP4, IRS1, AKT1, LEP, IGF1, IGF1R, EGF, GCG, RPS6KB1, INSR, PTPN1, IRS2, AKT2
GRP78 (HSPA5)	PDIA3, CALR, HSPA5, P4HB, HYOU1, DNAJC3, ATF6, CANX, PDIA4, HSP90B1, PDIA6, DNAJB11, HSPA1A, DNAJC10, DDIT3, ERN1, DNAJB1, XBP1, ATF4, PPIB
LY6E	LYPD1, LY6H, LYPD4, CEACAM5, LY6G6C, GPIHBP1, LYPD2, LY6G6D, TEX101, PSCA, LYPD5, LYPD3, LYPD6B, THY1, FOLR4, LY6E, MSLN, SPACA4, LYPD8, LY6K, LY6D
MADP1 (ZCRB1)	SF3B4, SF3B3, SART1, PRPF3, SNRNP200, SF3B2, PRPF31, SRSF7, SF3B1, HNRNPA1, SNRPD2, HNRNPA2B1, SF3B5, RNPC3, YBX1, ZCRB1, SNRNP35, PDCD7, SNRNP48, ZRSR2, ZMAT5, CDC5L, SF3A1, SNRPD3, PHF5A, SF3A2, SF3B6, SNRPB2, SNRPA1, SRSF1, PRPF6
TMPRSS2	FOXA1, CCND1, AKT1, TMPRSS2, AR, SRC, PTEN, CTNNB1, KLK3, TP53, PDGFRB, PIK3CA, HSP90AA1, PTK2, PIK3R1, NR3C1

Investigation of NAFLD-related genes and interacting genes/proteins of each SCEFs identified some shared genes (here determined as NS genes) (Table 2). There was no common gene between LY6E interacting genes/protein and NAFLD-related genes. Whereas, ACE2, ADAM17, DPP4, and TMPRSS2 had ≥ 9 interacting genes/protein existed among NAFLD-related genes. These 4 genes and their interacting genes/proteins were selected as hub genes.

Table 2

SCEFs and their interacting genes/proteins associated with NAFLD

SCEFs	Interacting genes associated with NAFLD
ACE2*	ACE, DPP4, GCG, REN, AGTR1, AGT, ACE2, APLN, CMA1
ADAM17*	TNF, IKBKG, IL10, EGFR, IKBKB, MAP3K7, TNFRSF1A, RIPK1, NOTCH1, ADAM17, JAG1, EGF, BIRC2
DPP4*	LEP, INS, INSR, AKT2, DPP4, GCG, IGF1, IRS2, AKT1, IRS1, AGT, PTPN1, EGF, RPS6KB1
TMPRSS2*	PTEN, FOXA1, PIK3CA, TP53, AKT1, CTNNB1, NR3C1, CCND1, HSP90AA1, AR, KLK3, PDGFRB
CTSL	CTSS, HLA-DQA1, HLA-DQB1, HLA-DQB2, HLA-DRB1
BSG	EGFR, ITGB1
CLEC4G	SYK, TREM1
CLEC4M	ROCK1
CTSB	CTSB, CTSD, CTSS, HLA-DQA1, HLA-DQB1, HLA-DQB2, HLA-DRB1
)	XBP1, ATF4, DDIT3, HSPA5, ERN1, PDIA3, HSPA1A
LY6E	−
MADP1 (ZCRB1)	HNRNPA1, HNRNPA2B1, PRPF31, SART1

*Top four SCEFs and their interacting genes were selected as hub genes

GO term enrichment analysis

The identified hub genes were investigated using Enrichr for GO and KEGG pathway analyses. The results of the GO analysis revealed that hub genes were involved in different proliferation and metabolism-related process, including regulation of cell proliferation, positive regulation of intracellular signal transduction, positive regulation of protein phosphorylation, positive regulation of cellular biosynthetic process, regulation of glucose import, cytokine-mediated signaling pathway, and regulation of glycogen biosynthetic biological processes (Fig. 1a). Moreover, KEGG pathway analysis revealed that the hub genes were highly associated with various pathways including cancers, insulin resistance, FoxO signaling pathway, human papillomavirus infection, Hepatitis C, PI3K-Akt, AMPK, and mTOR signaling pathways (Fig. 1b).

Fig. 1

Gene set enrichment analysis. Biological process (a) and KEGG pathway enrichment analysis of candidate hub genes (b)

Upstream regulatory network

Using the X2K online database, the TFs, PKs, and intermediate proteins regulating the hub genes were identified (Table 3). Also, the upstream regulatory network of the hub genes was constructed between the top ten TFs, top ten PKs, and their respective intermediate proteins (Fig. 2). Our analysis showed that the TFs and PKs were linked to 36 intermediate proteins. We identified AKT1 (connection > = 21) and MAPK14 (connection 10–14) as the most important PKs regulating hub genes expression. In addition, we classified TP53, AR, SP1, RELA, NFKB1, as the important TFs, since they shared the highest number of connections.

Table 3

Top ten TFs and PKs regulating the hub genes

	Adj P-value	Genes
TF
HIF1A	2.20E-11	ACE2, AR, ACE, CCND1, PIK3CA, LEP, AGTR1, KLK3, NR3C1
SP1	4.68E-11	IL10, EGF, INSR, PTEN, IRS2, KLK3, TNF, EGFR, AGT, AR, ADAM17, CCND1, LEP, AGTR1
RELA	5.11E-11	IL10, AR, CCND1, AKT2, PTEN, AKT1, KLK3, TNF, TP53, EGFR, AGT, BIRC2
NFKB1	5.11E-11	IL10, AR, CCND1, AKT2, PTEN, AKT1, KLK3, TNF, TP53, EGFR, AGT, BIRC2
EGR1	4.17E-10	AR, ACE, CCND1, PTEN, KLK3, TNF, TP53, EGFR
STAT3	1.72E-08	IL10, AR, CCND1, LEP, PTEN, AKT1, TP53, EGFR
AR	2.60E-08	FOXA1, PTPN1, PTEN, AKT1, CTNNB1, KLK3, EGFR
PGR	4.31E-08	IL10, CCND1, IRS2, KLK3, EGFR
PPARG	1.05E-07	JAG1, CCND1, PTEN, REN, TP53, EGFR
TP53	8.43E-07	PDGFRB, CCND1, PTEN, AKT1, CTNNB1, TP53, EGFR
PK
CSNK2A2	9.21E-08	IKBKB, PTPN1, HSP90AA1, IRS1, PTEN, CTNNB1, TNF, TP53, EGFR
ATM	9.21E-08	IKBKB, HSP90AA1, PTEN, AKT1, CTNNB1, IKBKG, NR3C1, TP53
CSNK2A1	2.30E-07	PTPN1, HSP90AA1, ACE, IRS1, PTEN, AKT1, CTNNB1, TNF, TP53
MAPK14	2.30E-07	IKBKB, PTPN1, AR, ADAM17, IRS1, AKT1, IRS2, NR3C1, TP53, EGFR
AKT1	2.99E-07	PTPN1, AR, IRS1, PTEN, AKT1, CTNNB1, IRS2, NR3C1
PRKCZ	1.40E-06	IKBKB, NOTCH1, IRS1, PTEN, AKT1
PRKDC	1.67E-06	IKBKB, HSP90AA1, IRS1, AKT1, IRS2, NR3C1, TP53
CHUK	2.40E-06	IKBKB, CCND1, IRS1, CTNNB1
MAPK3	2.85E-06	AR, ADAM17, IRS1, AKT1, NR3C1, TP53, EGFR
PRKCA	3.60E-06	IKBKB, PTPN1, AR, ADAM17, IRS1, INSR, AKT1, TP53, EGFR

Fig. 2

Upstream gene regulatory network of hub genes. Top ten TFs, top ten PKs, and related intermediate proteins are shown

miRNA-target gene interactions

miRTarBase_2017 was used to find the miRNAs that targeted the candidate hub gens. miRNA analysis discovered several miRNAs that potentially targeted the hub genes. 10 miRNAs with the most interactions to the target genes were selected as the top ten miRNAs (Table 4). Among them, hsa-miR-34a-5p contributed to the regulation of the highest number of target genes (n = 13). hsa-miR-34a-5p, hsa-miR-19a-3p and hsa-miR-152-3p regulated 13, 7 and 6 genes, respectively. CCND1 was the key gene that was targeted with 6 miRNAs. JAG1, PTEN, AKT1 were among the genes targeted by 4 miRNAs (Fig. 3).

Table 4

Top ten miRNAs regulating the expression of hub genes

miRNA	Adj P-value	Targets
hsa-miR-34a-5p	2.06E-06	IL10, PDGFRB, NOTCH1, JAG1, TNF, TNFRSF1A, CCND1, AGTR1, AKT1, CTNNB1, MAP3K7, TP53, BIRC2
hsa-miR-152-3p	4.93E-04	HSP90AA1, ADAM17, CCND1, IRS1, PTEN, APLN
hsa-miR-708-5p	0.017992059	CCND1, AKT2, AKT1, IKBKG
hsa-miR-200c-3p	0.017992059	IKBKB, NOTCH1, JAG1, IRS1, PTEN
hsa-miR-875-5p	0.017992059	CTNNB1, MAP3K7, EGFR
hsa-miR-143-3p	0.017992059	JAG1, AKT2, AKT1, TNF, MAP3K7
hsa-miR-145-5p	0.017992059	ADAM17, JAG1, IRS1, IRS2, EGFR
hsa-miR-19a-3p	0.017992059	IL10, CCND1, PIK3CA, PTEN, AKT1, TNF, TP53
hsa-miR-224-5p	0.017992059	PDGFRB, HSP90AA1, CCND1, APLN
hsa-miR-29a-3p	0.020130639	PDGFRB, CCND1, AKT2, PTEN, IGF1

Fig. 3

The regulatory network of top ten miRNA and their target genes

Discussion

While the world is suffering from the COVID-19 pandemic, there are groups of patients with comorbidities who are at an increased risk of hospitalization and even more susceptible to severe stages of infection. Among these comorbidities, NAFLD has shown to be one of the most important causes which significantly increases the odds of admission. A recent study on 6700 adults in the US proved that NAFLD/NASH is considerably one of the comorbidities with more risk of admission than age, gender, obesity, which raises the chances of hospitalization. It also showed patients who previously received metabolic syndrome treatment such as metformin had a notable drop in hospital admission rate after getting infected by SARS-CoV-2 (Bramante et al. 2020).

Here, we studied the NAFLD-related genes and SCEFs to find a link between molecular factors involved in NAFLD and COVID-19. Our analysis demonstrated an association of NAFLD-related genes with the most important SCEFs, including ACE2, ADAM17, TMPRSS2, and DPP4, which might describe the correlation for the higher severity ratio of SARS-CoV-2 infected patients with pre-existing NAFLD.

Construction of PPI network for each SCEF, MCODE and Cytohubba analysis were used to select more important genes interacting with 12 selected SCEFs. We found that each SCEF interacted with different genes/proteins. The common genes between NAFLD-related genes and SCEF-related genes were considered as the hub genes, which were more analyzed with different bioinformatic approaches to find different molecular mechanisms controlling the hub genes expression. The results of GO analysis of hub genes illustrated that most of them are enriched in biological processes such as regulation of cell proliferation, positive regulation of intracellular signaling transduction, and positive regulation of protein phosphorylation. GO analysis shows that AKT1 (as an interactor associated with NAFLD) correlates with viral entry factors TMPRSS2 and DPP4. Interestingly, AKT1 is involved in some cell proliferation and metabolism controlling processes such as regulation of cell proliferation, positive regulation of protein phosphorylation, regulation of glycogen biosynthetic, glycogen metabolic processes, and glucose import. It has been reported that AKT1 exhibits a great impact on viral gene transcription and little effect on viral protein synthesis in cells overexpressing the constitutively active AKT1 (Wang et al. 2014). Moreover, targeting AKT and PI3K/Akt/mTOR pathway as a treatment in COVID-19 patients has been suggested to protect lung epithelial cells and reduce fibrogenesis (Xia et al. 2020). AKT1 and AKT2 are also involved in Insulin resistance and FoxO signaling pathways. Dysregulated PI3K/AKT pathway in hepatocytes is commonly present in NAFLD (Matsuda et al. 2013). Based on these evidences, we can suggest that a high level of AKT1 expression may be involved in the severity of COVID-19 symptoms in NAFLD patients.

Here we found that NAFLD-related genes are also involved in cell entry mechanisms of SARS-CoV-2. Among 12 selected SCEFs, ACE2, ADAM17, TMPRSS2, and DPP4 had the highest number of interacting factors participated in NAFLD. ACE2 and TMPRSS2, known to be the co-receptors of SARS-CoV-2, are expressed in multiple organs and different tissues which may assist viral entry into host cells (South et al. 2020). Surface expression of ACE2 can affect the pathogenesis of COVID-19 and its blocking has the potential to stop viral entry and replication (Palau et al. 2020). In a recent study, overexpression of ACE2 has been shown to facilitate SARS-CoV-2 replication and greatly impact the infection in elderly and hypertensive patients (Peron and Nakaya 2020). TMPRSS2 is an androgen-responsive serine protease, expressed in multiple tissues and proved to participate in SARS-CoV-2 entry (Strope et al. 2020). TMPRSS2 and its related proteases can cleave ACE2 and SARS-CoV S protein (SARS-S). TMPRSS2-mediated ACE2 cytoplasmic tail cleavage may increase viral uptake and cleaved SARS-S may activate S protein used in membrane fusion (Heurich et al. 2014). On the other hand, multi-organ expression of ACE2 and TMPRSS2 has been seen in NAFLD patients and increased hepatic expression of ACE2 and TMPRSS2 in individuals with NAFLD may be a cause of hepatic complications in COVID-19 patients with NAFLD background (Lu et al. 2020; Singh et al. 2020; Meijnikman et al. 2021).

ADAM17 has been proved to facilitate entry of SARS-CoV-2 and induce tissue damage through TNF-α (Palau et al. 2020; Zipeto et al. 2020). It Is associated with JAG1 as an interactor associated with NAFLD. An overexpression in JAG1 has been seen in advanced NAFLD patients (Hotta et al. 2017). ADAM17 is involved in the regulation of cell proliferation and cytokine-mediated signaling pathways. A recent study showed that inhibition of ADAM17 may have a protective effect on SARS-CoV-2 infection by reducing the viral load (Palau et al. 2020). Moreover, tumor necrosis factor (TNF) and IL-10 that are important pro-inflammatory cytokines interact with ADAM17. TNF influences cytokine-mediated signaling, mTOR signaling, and insulin resistance pathways. Also, IL-10 is involved in the regulation of cell proliferation, cytokine-mediated signaling, and FoxO signaling pathways. In a meta-analysis including 50 studies with 7,865 patients, it is proved that there is a significant rise in TNF-α and IL-10 in the severe group compared to the non-severe group of COVID-19 patients (Akbari et al. 2020) and high levels of IL-10 in COVID-19 patients may lead to exhaustion of T-cells which limits their function (Diao et al. 2020). It has been reported that anti-TNF therapy shows a better prognosis in COVID-19 infection. Expression of TNF in NAFLD may also be an underlying reason for the severity of COVID-19 (Robinson et al. 2020).

DPP4 is another important SCEF that was identified in our study. It is a glycoprotein of 110 kDa, which is ubiquitously expressed on the surface of a variety of cells. It has been shown that high levels of plasma DPP4 are a predictor of the onset of metabolic syndrome, and DPP4 upregulation may be a possible underlying cause of COVID-19 disease severity (Bassendine et al. 2020). DPP4 is also expressed in several tissues, including the respiratory tract which facilitates SARS-CoV-2 entry, and by inducing cytokine storm may lead to fatal pneumonia (Solerte et al. 2020).

It has been shown that COVID-19 patients facing a severe form of the disease have experienced a cytokine storm occurrence. Among the top ten identified TFs, NFKB1 and STAT3 have been shown to take part in this hyper-inflammatory state. These TFs can regulate proteins such as IL-6 and TNF that cause a cytokine storm in the host cells (Hojyo et al. 2020). We also identified the top ten PKs. Among them, MAPK14 has been reported to be a potential target in COVID-19 patients. In a previous study to explore the mechanism of therapeutic effect of vitamin A in the treatment of COVID-19, It has been reported that MAPK14 is one of the core targets of vitamin A against COVID-19 (Li et al. 2020).

In early 2020, scientists and physicians started studying different antiviral and anticancer drugs that might decrease COVID-19 pandemic death rates. Remdesivir is a broad-spectrum antiviral drug that has shown to be a therapeutic option for COVID-19 by resembling an adenosine triphosphate (ATP) molecule and competing with the nucleotide in the synthesis of the viral RNA (Saha et al. 2020). Interestingly, there is evidence that prove Remdesivir reduces the inflammatory response in high-fat diet-induced NAFLD in mice that had notably increased immune cell infiltration and inflammatory response. These inflammatory effects were triggered by high-fat diet-induced NAFLD and showed an increase in the level of several cytokines such as NFkB, TBK1, IRF3, IFN-b, TNF-a, IL-6, and IL-18 (Li and Su 2020). It has also been reported that Andrographolide, a potential inhibitor of SARS-CoV-2, can bind to NFkB1 and block TNF-induced cytokine storm in COVID-19 (Rehan et al. 2021).

miRNAs function in RNA silencing and post-transcriptional regulation of gene expression. Here, we identified the top 10 miRNAs that target the selected hub genes. Among these miRNAs, miR-34a, miR-29a, miR-200c, and miR-152 were proven to have been implicated in liver diseases such as NAFLD, by playing their regulatory roles in liver metabolism (Feng et al. 2014; Gjorgjieva et al. 2019). Moreover, it has been affirmed that SARS-CoV-2 disturbs the host immune system by manipulating its miRNAs. Additionally, this virus encodes its miRNAs and deactivate gene function in the host cells (Zhang et al. 2021). A recent study has shown that miR-19a-3p was upregulated in the plasma of COVID-19 patients, indicating its role in the SARS-CoV-2 infection (Fayyad-Kazan et al. 2021). It has been demonstrated that measuring the expression levels of miR-29a can provide essential predictive (López-Riera et al. 2017; Jampoka et al. 2018; Lambrecht et al. 2019) and diagnostic (Jampoka et al. 2018) information about NAFLD, as well as its progression (Lin et al. 2019). This microRNA mediates biological functions, including epigenetic modification, neutralizing oxidative stress, exerting an anti-inflammatory effect on the pathogenesis of NAFLD, regulating mitochondrial metabolism, and preventing lipid accumulation in the liver, which can protect the liver from developing damage and consequently dysfunction (Lin et al. 2020). A recent study has shown that miR-29a which is found significantly in higher amounts in the blood of COVID-19 patients, is suggested as a potential biomarker for diagnosis and monitoring the disease (Donyavi et al. 2021). miR-29a is a predicted miRNA binder to SARS-CoV-2 genome in two regions of ORF1ab and ORF9 which is highly expressed in the lungs of COVID-19 patients (Pierce et al. 2020). According to our data, hsa-miR-34a-5p is an important gene regulator among the miRNAs discovered here, that targets 13 genes including CCND1, AKT1, MAP3K7, IL10, PDGFRB, NOTCH1, JAG1, TNF, TNFRSF1A, AGTR1, CTNNB1, TP53, and BIRC2. When miR-34a is overexpressed, it increased steatosis as it directly targets downregulation of hepatic PPARα and SIRT1; and interestingly, when it’s inhibited, AMPK pathway is activated and lipid accumulation is suppressed (Castro et al. 2013; Ding et al. 2015; Wang et al. 2020b). Circulating miR-34a found in higher serum levels of NAFLD patients can be a biomarker of the disease, as well (Yamada et al. 2013). Moving on to COVID-19 findings in correlation with miR-34a, post-mortem lung biopsies of COVID-19 patients depicted reduced expression of miR-34a which based on functional enrichment analysis is associated with endothelial dysfunction and inflammatory response (Centa et al. 2020).

Conclusions

To conclude, one of the most vulnerable groups of the population feared to manifest severe symptoms of SARS-CoV-2 infection are amongst those with underlying health conditions including NAFLD patients. In this research, we investigated several genetic factors that could determine the reason behind the higher disease severity of COVID-19 in patients with NAFLD background. We identified AKT1, IL10, TNF, and several other key proteins in NAFLD patients as potential risk factors. Moreover, regulation of several PKs and TFs such as MAPK14, HIF1A, SP1, and NFKB1, along with cytokine-mediated signaling, PI3K-Akt signaling pathway, and miRNAs enriched in NAFLD were also found to affect the progression and severity of COVID-19. We hope our study sheds light on developing COVID-19 therapeutics and discovering the causes of the severity of SARS-CoV-2 infection in NAFLD patients.