Deciphering the link between Diabetes mellitus and SARS-CoV-2 infection through differential targeting of microRNAs in the human pancreas.

PURPOSE: Coronavirus Disease 2019 (COVID-19) severity and Diabetes mellitus affect each other bidirectionally. However, the cause of severe acute respiratory syndrome coronavirus 2 (SARS-CoV2) infection on the incidence of diabetes is unclear. In the SARS-CoV-2-infected cells, host microRNAs (miRNAs) may target the native gene transcripts as well as the viral genomic and subgenomic RNAs. Here, we investigated the role of miRNAs in linking Diabetes to SARS-CoV-2 infection in the human pancreas.
METHODS: Differential gene expression and disease enrichment analyses were performed on an RNA-Seq dataset of human embryonic stem cell-derived (hESC) mock-infected and SARS-CoV-2-infected pancreatic organoids to obtain the dysregulated Diabetes-associated genes. The miRNA target prediction for the Diabetes-associated gene transcripts and the SARS-CoV-2 RNAs has been made to determine the common miRNAs targeting them. Minimum Free Energy (MFE) analysis was done to identify the miRNAs, preferably targeting SARS-CoV-2 RNAs over the Diabetes-associated gene transcripts.
RESULTS: The gene expression and disease enrichment analyses of the RNA-Seq data have revealed five biomarker genes, i.e., CP, SOCS3, AGT, PSMB8 and CFB that are associated with Diabetes and get significantly upregulated in the pancreas following SARS-CoV-2-infection. Four miRNAs, i.e., hsa-miR-298, hsa-miR-3925-5p, hsa-miR-4691-3p and hsa-miR-5196-5p, showed preferential targeting of the SARS-CoV-2 genome over the cell's Diabetes-associated messenger RNAs (mRNAs) in the human pancreas.
CONCLUSION: Our study proposes that the differential targeting of the Diabetes-associated host genes by the miRNAs may lead to diabetic complications or new-onset Diabetes that can worsen the condition of COVID-19 patients.

Introduction

Coronavirus Disease (COVID-19) is a fast-spreading disease that has caused a global crisis. This pandemic is a highly infectious viral disease caused by Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [1, 2]. According to the World Health Organization’s (WHO’s) COVID-19 Dashboard [3], 200,840,180 confirmed COVID-19 cases resulting in 4,265,903 deaths have been reported globally as of 6th August 2021.

Diabetes is the chief cause of death in the hospitalized COVID-19 patients [4]. Patients having Diabetes mellitus (commonly referred to as Diabetes) of either Type 1 or Type 2 show poor prognosis with SARS-CoV-2 infection due to blood glucose level fluctuations and metabolic complications. An increase in blood glucose level helps in the SARS-CoV-2 replication and proliferation in the human monocytes [5, 6]. New-onset diabetes has also been observed after COVID-19 infection [7]. Thus, COVID-19 severity and Diabetes affect each other in a bidirectional manner.

Moreover, diabetic patients are also at increased risk of COVID-19. The resultant hyperglycemia can diminish the individual’s immune response towards handling the viral infection. COVID-19 mortality is also amplified for diabetic patients due to complications like cardiovascular diseases and kidney-related problems [7-10]. Lately, it has been reported that poorly controlled Diabetes mellitus in COVID-19 patients may lead to Mucormycosis, which is becoming an area of rising concern [11, 12]. Thus, regulating the blood glucose levels and preventing diabetic complications is required for diabetic COVID-19 patients to avert severe consequences of the infection [5]. Autopsy reports and experimental evidence have been known to reveal SARS-CoV-2 infection in the human pancreas, including the insulin-producing islet beta cells [13-16]. It is possible due to the angiotensin-converting enzyme 2 (ACE2) receptors on the pancreas cells. The SARS-CoV-2 virus can reach the pancreas from the duodenal epithelium [17]. ACE2 acts as a high-affinity receptor on the cell for the spike protein of SARS-CoV-2 [17-19]. The expression of ACE2 and other SARS-CoV-2 entry factors, i.e., Transmembrane serine protease 2 (TMPRSS2), Neuropilin 1 (NRP1) and Transferrin receptor protein 1 (TRFC) in the insulin-producing beta cells of the pancreas pave the way for the infection in these cells [14, 16].

SARS-CoV-2 has a positive-sense single-stranded ribonucleic acid (+ ssRNA) genome enclosed in a protein-containing lipid bilayer [5, 20]. In the host cell, the SARS-CoV-2 genome is translated to produce polypeptides which are further cleaved to form non-structural proteins. The non-structural proteins help in the replication and transcription of the viral genome. The positive-sense genomic RNA and the subgenomic RNAs are synthesized from the negative-sense RNA intermediate template. The subgenomic RNAs are translated to form the structural and accessory proteins of the virus. The positive-sense viral genome copies are packaged to form virions by the structural proteins, i.e., spike protein [S], envelope protein [E], membrane protein [M], and nucleocapsid protein [N]. The accessory proteins, i.e., 3a, 6, 7a, 7b, 8, and 10, are translated from their respective open reading frames with the same names. The subgenomic RNAs share a common leader five prime untranslated regions (5ʹUTR) and three prime untranslated regions (3ʹUTR) sequences with the viral genome RNA [21, 22]. In the host cell, the native messenger RNAs (mRNAs) translation is more efficient than the SARS-CoV-2 genome, but the high number of viral transcripts dominates virus translation [22].

MicroRNAs (miRNAs) play an essential role in the regulation of gene expression in a cell. In a virus-infected cell, the host cell’s miRNAs help in cell defense by targeting the viral RNA genome and its transcribed RNAs. In humans, the miRNAs primarily target the 3ʹUTR of the mRNAs. However, the miRNAs primarily target the 3ʹUTR and 5ʹUTR of the viral RNA genomes. Thus, during a viral infection, the host cell’s miRNAs may target the viral genome rather than the host mRNAs. It leads to a competition for miRNA regulation between the host cell’s mRNAs and the viral genome copies in the cell. This differential targeting of the miRNAs may result in the dysregulation of the host cell’s genes [23-26]. To decipher the complex relationship between Diabetes and COVID-19, we tried to understand the interplay between the SARS-CoV-2 RNAs, cell’s host mRNAs, and miRNAs in the infected pancreatic cells. For this purpose, we have analyzed the gene expression data of the mock-infected and SARS-CoV-2-infected human embryonic stem cell (hESC)-derived pancreatic organoids in this study.

Materials and methods

Differential gene expression analysis

For the study, the gene expression data were retrieved from the Gene Expression Omnibus (GEO), National Center for Biotechnology Information [27]. The GEO dataset accession number, GSE151803, includes the RNA-Seq data of the hESC liver, lung, and pancreatic cells/organoids obtained using the Illumina NovaSeq 6000 sequencing system. The dataset has been designed for transcriptomic analysis of the mock-infected and SARS-CoV-2-infected organoids [28]. The RNA-Seq data of three mock-infected and three SARS-CoV-2-infected hESC pancreatic organoids present in the dataset were selected for this study. For obtaining the pancreatic organoids, hESC line-RUES2 was cultured, followed by performing pancreatic endocrine cell differentiation in the tissue. The Differentially Expressed Genes (DEGs) between the mock-infected and the SARS-CoV-2-infected pancreatic organoids were identified using the DESeq2 R package [29, 30]. The genes with the adjusted p value < 0.05 and the absolute log2 fold change value > 1 (padj < 0.05 and |log2FC|> 1) were considered to be the significant DEGs. Furthermore, single-cell RNAseq (scRNA-seq) data of COVID-19-infected human pancreatic organoid from gene expression omnibus datset GSE159556 was utilized to cross-validate our observation [31]. The analysis of scRNA-seq was performed using Seurat 4.0.2 package [32].Volcano plots were generated using Enhanced Volcano R package [33].

Disease enrichment analysis

Disease enrichment analysis based on the DEGs was done using DAVID-Functional Annotation Tool [34]. For this, the overrepresented disease terms associated with the DEGs were enriched from the Gene-Disease Associations Dataset (GAD) [35] linked to the tool. The resulting pancreas-associated disease terms and the DEGs associated with them were selected for further study. In addition, the disease-associated genes encoding for secretory proteins were checked using The Human Protein Atlas [36-38].

Gene–Gene Interaction network analysis

The Gene–Gene Interaction (GGI) network between the DEGs was constructed with the Cytoscape-GeneMANIA app [39-41]. The GGI network was visualized to decipher the associations between the DEGs like co-expression, genetic interactions, pathways, co-localization, and protein domain similarity. In addition, network topological analysis was also done to get an insight into the influence of hub DEG nodes on the network. Three network parameters, i.e., degree, closeness centrality and betweenness centrality of the gene nodes, were calculated using the Cytoscape-Network Analyzer plugin [42].

miRNA target prediction

The human miRNAs targeting the UTRs of the SARS-CoV-2 genome, hereafter referred to as the CoV-tar-miRNAs, were predicted using the miRDB online tool [43]. The complete genome reference sequence of SARS-CoV-2, Wuhan-Hu-1, was retrieved from NCBI RefSeq [44] ID, NC_045512.2 [45]. The 3ʹUTR and 5ʹUTR of the viral RNA genome were selected for the study. The target genes of CoV-tar-miRNAs were obtained using the Predicted Target Module of miRWalk2.0 [46, 47]. In the module’s input parameters, the minimum seed length was chosen to be seven and the p value cutoff was kept to be 0.05. The putative target genes of the CoV-tar-miRNAs have been obtained by selecting four databases for the results, i.e., miRWalk, RNA22, miRanda and Targetscan, all included in miRWalk2.0 [46, 47]. As a result, the Diabetes-associated DEGs in the target gene list of CoV-tar-miRNAs were obtained. Also, the common miRNAs between the CoV-tar-miRNAs and miRNAs targeting the Diabetes-associated DEGs were considered to be the Diabetes-associated CoV-tar-miRNAs. Furthermore, for verifying the availability of miRNAs in the normal human pancreas tissue, their expression values were checked in the TissueAtlas database [48] and GeneAnalytics tool of GeneCards Suite [49, 50].

Minimum free energy analysis

For evaluating the miRNA binding with the 3ʹUTR sequence of the target Diabetes-associated genes’ transcripts, the state-of-the-art prediction tool RNAhybrid was used. RNAhybrid predicts secondary structures between the miRNA and the target mRNA through Minimum Free Energy (MFE) calculations. [51, 52]. The 3ʹUTR nucleotide sequences of the Diabetes-associated genes’ transcripts were retrieved from NCBI.

Results

Differential gene expression and disease enrichment analysis

We identified 30 DEGs between the mock-infected and SARS-CoV-2-infected human pancreatic organoids, among which 26 are upregulated, and 4 are downregulated in the SARS-CoV-2-infected hESC pancreatic organoids (Table 1). The statistically significant DEGs, i.e., DEGs with padj < 0.05 and |log2FC|> 1, are shown with red dots in the presented volcano plot (Fig. 1a). The gene-based disease enrichment of the 30 DEGs resulted in “Type 1 Diabetes” as the most significant disease term. Four upregulated DEGs, i.e., CP, SOCS3, AGT, and PSMB8, were linked with Type 1 Diabetes. Also, two upregulated DEGs, i.e., CP and CFB, were enriched for the term “insulin” (Table 2). COVID-19 is linked with other enriched disease terms, but only these two terms were selected for the study due to their direct association with the pancreas. The Human Protein Atlas [30-32] revealed that among the DEGs associated with enriched disease terms, CP, AGT, CFB, SERPINA3, CXCL2, C8B, and AKR1B10 encode for secretory proteins. We cross-validated our results by analyzing the single-cell RNAseq (scRNA-seq) data of COVID-19 infected human pancreatic organoid from gene expression omnibus dataset GSE159556 [31]. We found upregulation of AGT, CFB, PSMB8 in acinar cells and ductal cells, however, SOCS3 was upregulated in beta cells (Fig. S1.A, B, C). We also checked the expression of AGT, CFB, PSMB8, SOCS3, and CP in different pancreas cell types, namely alpha cells, acinar cells, beta cells, delta cells, ductal cells, pp cells, endothelial cells, mesenchyme cells (Fig. S1 D).

Table 1

DEGs between mock-infected and SARS-CoV-2-infected hESC pancreatic organoids

S. No	Gene symbol	Gene description	Log2FoldChange
1	CP	Ceruloplasmin	3.218646021
2	SOCS3	Suppressor of cytokine signaling 3	2.897737916
3	VNN3	Vanin 3	2.676974204
4	CEBPD	CCAAT enhancer binding protein delta	2.041519486
5	FAM169B	Family with sequence similarity 169 member B	1.887720767
6	VNN2	Vanin 2	1.650151224
7	HPX	Hemopexin	1.576405664
8	EVA1A	Eva-1 homolog A, regulator of programmed cell death	1.483037767
9	SERTM1	Serine rich and transmembrane domain containing 1	1.482242018
10	CXCL2	C-X-C motif chemokine ligand 2	1.464664038
11	C8B	Complement C8 beta chain	1.45102627
12	CFB	Complement factor B	1.449773368
13	SERPINA3	Serpin family A member 3	1.282619509
14	TACR1	Tachykinin receptor 1	1.253706268
15	HCP5	HLA complex P5	1.228387286
16	CHRD	Chordin	1.226960717
17	GADD45G	Growth arrest and DNA damage inducible gamma	1.223936929
18	APOL6	Apolipoprotein L6	1.164596924
19	SAMD11	Sterile alpha motif domain containing 11	1.129551911
20	AGT	Angiotensinogen	1.106392771
21	C9orf16	Chromosome 9 open reading frame 16	1.096117015
22	NLRC5	NLR family CARD domain containing 5	1.094272111
23	SUCNR1	Succinate receptor 1	1.067266654
24	PSMB8	Proteasome 20S subunit beta 8	1.031089013
25	UNC5CL	Unc-5 family C-terminal like	1.015471709
26	MAP6D1	MAP6 domain containing 1	1.000047502
27	AKR1B10a	Aldo–Keto Reductase Family 1 member B10	− 1.942847572
28	SYNDIG1La	Synapse differentiation inducing 1 like	− 1.628255609
29	PKHD1L1a	Polycystic kidney and hepatic disease 1 like 1	− 1.414217069
30	TMEM236a	Transmembrane protein 236	− 1.040226661

aGenes downregulated in SARS-CoV-2-infected hESC pancreatic organoids

Fig. 1

a Volcano plot showing the significant DEGs between mock-infected and SARS-CoV-2-infected hESC pancreatic organoids. Differentially expressed genes with padj < 0.05 and |log2FC|> 1 are shown in red dots b GGI network between the DEGs: Diabetes-associated (Type 1 Diabetes and insulin-associated) DEGs are connected directly or indirectly through co-expression c Diabetes-associated gene targets of the miRNAs that can target the SARS-CoV-2 genome and the Diabetes-associated gene: The circular nodes depict the Diabetes-associated genes. The elliptical nodes represent the miRNAs that can target the 3ʹUTR of the SARS-CoV-2 genome and the Diabetes-associated genes. The rectangular nodes represent the miRNAs that can target the 5ʹUTR of the SARS-CoV-2 genome and the Diabetes-associated genes

Table 2

Gene-based disease enrichment of DEGs

S. No	Disease term	Genes
1	Type 1 diabetes	SOCS3, AGTa, CPa, PSMB8
2	Alzheimer’s disease, Parkinson’s disease, insulin, lung function, depression, longevity	CFBa, CPa
3	Cerebrovascular disease	AGTa, SERPINA3a
4	Respiratory syncytial virus bronchiolitis, asthma, bronchiolitis	SOCS3, CXCL2a, PSMB8
5	Birth weight, leukemia, acute myeloid leukemia, precursor cell lymphoblastic leukemia-lymphoma, meningeal neoplasms, meningioma, non-hodgkin’s lymphoma	C8Ba, SOCS3, CFBa
6	Nephropathy	AGTa, AKR1B10a
7	Myocardial Infarction	EVA1A, AGTa, GADD45G, SERPINA3a
8	Plasma HDL cholesterol (HDL-C) levels	SOCS3, CEBPD, AGTa
9	Macular degeneration	C8Ba, CFBa, SERPINA3a

aGenes encoding for secretory proteins

GGI network analysis

The GGI network between the 29 DEGs was constructed to determine the types of relations between the DEGs. The GGI network revealed that the Diabetes-associated DEGs, i.e., CP, SOCS3, AGT, PSMB8, and CFB, are connected directly or indirectly through co-expression (Fig. 1b, Supplemental Table S1). The network excludes the HCP5 gene as it codes for a long non-coding RNA, and its data are not present in GeneMANIA [50]. In the network, all the Diabetes-associated DEGs show direct or indirect co-expression interactions between them. The CP-AGT, CP-CFB, SOCS3-CFB, and AGT-CFB node pairs have a direct co-expression interaction between them. CP-AGT, CP-CFB, and AGT-CFB node pairs also show co-localization interaction. CP-SOCS3 and SOCS3-AGT node pairs have an indirect co-expression interaction through CFB. PSMB8-CFB node pair indirectly interacts through CP through co-expression. CP-PSMB8 node pair shows indirect co-expression interaction between them through APOL6. SOCS3-PSMB8 and AGT-PSMB8 node pairs show indirect co-expression interaction through the CFB-CP-APOL6 and CP-APOL6 gene paths, respectively. Two downregulated DEGs, i.e., AKR1B10 and PKHD1L1 show direct interactions with the upregulated DEGs (Fig. 1b, Supplemental Table S2).

Furthermore, to decipher the influential genes in the DEGs’ network, three topological parameters, i.e., degree, closeness centrality, and betweenness centrality, were calculated for the gene nodes (Table 3). The Diabetes-associated gene CP with 11 links depicts the hub gene with the highest closeness and betweenness centralities, and thus, it is the most influential gene in the network. CFB gene is the second-highest influential and diabetes-associated gene in terms of degree (10 links) and closeness centrality values.

Table 3

Topological parameters of the gene nodes in GGI network (in descending order of Closeness Centrality)

S. No	Gene name	Degree	Closeness Centrality	Betweenness Centrality
1	SUCNR1	1	1	0
2	TACR1	1	1	0
3	CPb	11	0.65625	0.334650416
4	SERPINA3	9	0.6	0.078909675
5	CFBb	10	0.567567568	0.124459562
6	C8B	8	0.538461538	0.124801587
7	EVA1A	7	0.525	0.007420635
8	VNN3	6	0.525	0.055616024
9	VNN2	6	0.525	0.25037037
10	AGTb	7	0.512195122	0.062619048
11	CXCL2	8	0.512195122	0.101135676
12	HPX	6	0.5	0
13	CHRD	3	0.446808511	0.034950869
14	APOL6	2	0.4375	0.096216931
15	GADD45G	4	0.4375	0.022896825
16	PKHD1L1	1	0.403846154	0
17	SOCS3b	3	0.396226415	0.007539683
18	UNC5CL	2	0.381818182	0.002380952
19	FAM169B	2	0.375	0.07521164
20	CEBPD	3	0.368421053	0.002380952
21	AKR1B10	1	0.355932203	0
22	SAMD11	1	0.35	0
23	PSMB8b	2	0.328125	0.023994709
24	NLRC5	2	0.291666667	0.013492063
25	C9orf16a	0	0	0
26	TMEM236a	0	0	0
27	SYNDIG1La	0	0	0
28	SERTM1a	0	0	0
29	MAP6D1a	0	0	0

aSingle node in the network

bDifferentially expressed genes associated with diabetes

miRNAs targeting the 3ʹUTR and 5ʹUTR of the viral genome

The SARS-CoV-2 genome consists of a linear 29,903 nucleotides long ssRNA [45]. The 3ʹUTR of the genome is 229 nucleotides long, including a polyA tail. It lies in the viral genome at the position from 29,675 to 29,903 nucleotides. The 5ʹUTR of the genome is 265 nucleotides long and spans from 1 to 265 nucleotides in the viral genome (Supplemental Table S3). We conducted nucleotide sequence-based miRNA target prediction analysis on the SARS-CoV-2 genome UTRs to identify the CoV-tar-miRNAs. Eleven miRNAs, i.e., hsa-miR-298, hsa-miR-7851-3p, hsa-miR-1303, hsa-miR-3925-5p, hsa-miR-8075, hsa-miR-4691-3p, hsa-miR-1283, hsa-miR-3123, hsa-miR-5196-5p, hsa-miR-4747-5p and hsa-miR-4645-3p have been predicted as 5ʹUTR CoV-tar-miRNAs, potentially targeting the 5ʹUTR of the SARS-CoV-2 genome (Table 4, Supplemental Table S4). Similarly, ten miRNAs, i.e., hsa-miR-3941, hsa-miR-466, hsa-miR-4775, hsa-miR-4717-3p, hsa-miR-5088-5p, hsa-miR-603, hsa-miR-6749-3p, hsa-miR-1236-3p, hsa-miR-4279 and hsa-miR-4672 have been predicted as the 3ʹUTR CoV-tar-miRNAs that can potentially target the 3ʹUTR of the viral genome (Table 4, Supplemental Table S5).

Table 4

miRNAs potentially targeting the viral genome, i.e., CoV-tar-miRNAs

S. No	5ʹUTR CoV-tar-miRNAs	3ʹUTR CoV-tar-miRNAs
1	hsa-miR-298	hsa-miR-3941
2	hsa-miR-7851-3p	hsa-miR-466
3	hsa-miR-1303	hsa-miR-4775
4	hsa-miR-3925-5p	hsa-miR-4717-3p
5	hsa-miR-8075	hsa-miR-5088-5p
6	hsa-miR-4691-3p	hsa-miR-603
7	hsa-miR-1283	hsa-miR-6749-3p
8	hsa-miR-3123	hsa-miR-1236-3p
9	hsa-miR-5196-5p	hsa-miR-4279
10	hsa-miR-4747-5p	hsa-miR-4672
11	hsa-miR-4645-3p

Human miRNAs targeting the 3ʹ and 5ʹ UTR of the SARS-CoV-2 genome and Diabetes-associated DEGs

To find the CoV-tar-miRNAs potentially targeting the Diabetes-associated DEGs, their gene targets were identified. First, the set of Diabetes-associated DEGs and the miRNAs potentially targeting them were obtained from the list of gene targets. Then, the miRNAs common in potentially targeting the CoV-tar-miRNAs and Diabetes-associated DEGs were obtained. These Diabetes-associated CoV-tar-miRNAs can potentially target the SARS-CoV-2 UTRs as well as the transcripts of Diabetes-associated DEGs. We found that five Diabetes-associated DEGs, i.e., CP, PSMB8, SOCS3, AGT, and CFB, can be targeted by one or more CoV-tar-miRNAs (Fig. 1c, Table 5). Eight 3ʹUTR and nine 5ʹUTR CoV-tar-miRNAs potentially target the Diabetes-associated DEGs. The CoV-tar-miRNAs: hsa-miR-466 and hsa-miR-4775 target the 3ʹUTR while hsa-miR-3123, hsa-miR-4691-3p, hsa-miR-4747-5p and hsa-miR-5196-5p target the 5ʹUTR of the CP mRNA. The CoV-tar-miRNAs: hsa-miR-1236-3p, hsa-miR-3941, hsa-miR-4279, hsa-miR-4717-3p, hsa-miR-5088-5p and hsa-miR-6749-3p target the 3ʹUTR, whereas hsa-miR-1303, hsa-miR-3123, hsa-miR-3925-5p, hsa-miR-4691-3p and hsa-miR-7851-3p target the 5ʹUTR of the SOCS3 mRNA. The CoV-tar-miRNA hsa-miR-4775 targets the 3ʹUTR. However, hsa-miR-1303, hsa-miR-298, hsa-miR-4645-3p and hsa-miR-7851-3p target the 5ʹUTR of the AGT mRNA. The CoV-tar-miRNAs: hsa-miR-1236-3p, hsa-miR-4279 and hsa-miR-5088-5p target the 3ʹUTR of the PSMB8 mRNA while hsa-miR-298, hsa-miR-3123 and hsa-miR-3925-5p target its 5ʹUTR. hsa-miR-4775 targets the 3ʹUTR of the CFB gene. hsa-miR-3123 and hsa-miR-4775 are the most influential CoV-tar-miRNAs, each targeting three upregulated Diabetes-associated genes. Thus, four CoV-tar-miRNAs commonly target SOCS3 and PSMB8. One CoV-tar-miRNA targets CP, AGT, and CFB. One CoV-tar-miRNA targets CP, SOCS3, and PSMB8. One CoV-tar-miRNA targets CP and SOCS3. One CoV-tar-miRNA targets AGT and PSMB8. (Fig. S2, Supplemental Table S6). We checked the expression data of the Diabetes-associated CoV-tar-miRNAs in the human pancreas across the published experimental work. Except for hsa-miR-466, hsa-miR-3123, hsa-miR-6749-3p, and hsa-miR-7851-3p, all the Diabetes-associated CoV-tar-miRNAs were found to be expressed in the human pancreas [48, 50, 53] (Table 5, Supplemental Table S7).

Table 5

microRNAs potentially targeting the Diabetes-associated DEGs and SARS-CoV-2 genome

S. No	Diabetes-associated DEGs	SARS-CoV-2 genome’s region targeted by CoV-tar-miRNAs	CoV-tar-miRNAs targeting the DEGs
1	CP	3ʹUTR	hsa-miR-466a hsa-miR-4775b
		5ʹUTR	hsa-miR-3123a hsa-miR-4691-3pb hsa-miR-4747-5pb hsa-miR-5196-5pb
2	SOCS3	3ʹUTR	hsa-miR-1236-3pb hsa-miR-3941b hsa-miR-4279b hsa-miR-4717-3pb hsa-miR-5088-5pb hsa-miR-6749-3pa
		5ʹUTR	hsa-miR-1303b hsa-miR-3123a hsa-miR-3925-5pb hsa-miR-4691-3pb hsa-miR-7851-3pa
3	AGT	3ʹUTR	hsa-miR-4775b
		5ʹUTR	hsa-miR-1303b hsa-miR-298b hsa-miR-4645-3pb hsa-miR-7851-3pa
4	PSMB8	3ʹUTR	hsa-miR-1236-3pb hsa-miR-4279b hsa-miR-5088-5p
		5ʹUTR	hsa-miR-298b hsa-miR-3123a hsa-miR-3925-5pb
5	CFB	3ʹUTR	hsa-miR-4775b

aUnknown expression of miRNA in the human pancreas tissue

bExpressed in the human pancreas[42, 44, 47]

Furthermore, to evaluate whether the Diabetes-associated CoV-tar-miRNAs’ favorably hybridize with their viral and host Diabetes-associated RNA targets, we employed RNAhybrid 2.1.2 tool. The negative Minimum Free Energy (MFE) calculations suggested a favorable interaction of CoV-tar-miRNAs (hsa-miR-298, hsa-miR-3925-5p, hsa-miR-4691-3p, and hsa-miR-5196-5p) with the viral UTRs than with Diabetes-associated DEGs’ transcripts (Table 6). hsa-miR-4691-3p and hsa-miR-5196-5p show lower MFE with the viral 5ʹUTR than with the CP mRNA’s 3ʹUTR. hsa-miR-3925-5p show lower MFE with the viral 5ʹUTR than with the SOCS3 mRNA’s 3ʹUTR. hsa-miR-298 and hsa-miR-3925-5p show lower MFE with the viral 5ʹUTR than with the PSMB8 mRNA’s 3ʹUTR. We retrieved the 3ʹUTRs of the Diabetes-associated mRNAs from NCBI (Supplemental Table S3) for MFE analysis.

Table 6

Minimum Free Energy (MFE) of CoV-tar-miRNAs with the SARS-CoV-2 UTRs and Diabetes-associated host cell’s mRNA targets

3ʹUTR CoV-tar-miRNAs	MFE (kcal/mol) of miRNAs with 3ʹUTR of viral genome	Diabetes-associated host cell’s mRNA targets of CoV-tar-miRNAs	MFE (kcal/mol) of miRNAs with 3ʹUTR of host cell’s mRNA targets
hsa-miR-3941	− 19.2	SOCS3	− 24.1
hsa-miR-466	− 14.7	CP	− 16.2
hsa-miR-4775	− 15.7	CP	− 17.3
		AGT	− 19.0
		CFB	− 18.7
hsa-miR-4717-3p	− 26.3	SOCS3	− 32.0
hsa-miR-5088-5p	− 21.8	SOCS3	− 34.6
		PSMB8	− 26.0
hsa-miR-6749-3p	− 21.4	SOCS3	− 35.6
hsa-miR-1236-3p	− 24.5	SOCS3	− 34.7
		PSMB8	− 25.7
hsa-miR-4279	− 19.3	SOCS3	− 31.8
		PSMB8	− 23.1
5ʹUTR CoV-tar-miRNAs	MFE of miRNAs with 5ʹUTR of viral genome
hsa-miR-298a	− 24.8	AGT	− 28.3
		PSMB8a	− 24.7
hsa-miR-7851-3p	− 26.7	SOCS3	− 30.7
		AGT	− 29.0
hsa-miR-1303	− 21.0	SOCS3	− 29.9
		AGT	− 25.5
hsa-miR-3925-5pa	− 24.2	SOCS3a	− 22.9
		PSMB8a	− 22.4
hsa-miR-4691-3pa	− 25.5	CPa	− 20.7
		SOCS3	− 32.1
hsa-miR-3123	− 14.7	CP	− 15.8
		SOCS3	− 17.0
		PSMB8	− 15.4
hsa-miR-5196-5pa	− 29.6	CPa	− 25.6
hsa-miR-4747-5p	− 21.3	CP	− 24.2
hsa-miR-4645-3p	− 20.9	AGT	− 29.0

aCoV-tar-miRNAs showing higher MFE with Diabetes-associated mRNA targets than with SARS-CoV-2 genome UTR

Discussion

Hospitalized COVID-19 patients show the highest mortality with diabetes as comorbidity [4]. The SARS-CoV-2 infection has been detected in different human organs, including the pancreas [13, 54]. In this work, we have investigated the role of human pancreas miRNAs in linking Diabetes to COVID-19. We have noted that in the SARS-CoV-2-infected hESC pancreas tissue, 26 and 4 genes were upregulated and downregulated, respectively. Among them, we have found that four upregulated genes, i.e., CP, SOCS3, AGT and PSMB8 are associated with Type 1 Diabetes while two upregulated genes are associated with the term “insulin”. Since insulin is associated with both Type 1 and Type 2 Diabetes, the CP and CFB genes can be associated with both the types of Diabetes. Among the genes related to other enriched diseases, six upregulated genes, i.e., CP, AGT, CFB, SERPINA3, CXCL2, and C8B, encode proteins secreted to blood [36-38]. These proteins may also reach other organs of the body, thus, being involved in other diseases. Therefore, we suggest that their roles in other diseases must be further investigated. For example, Diabetes can result in Nephropathy or Macular Degeneration [55-57]; thus, they can be the indirect results of COVID-19.

SARS-CoV-2 is a + ssRNA genome virus, and its genome itself acts as an mRNA. In addition to its replicated RNA genome copies, subgenomic RNAs are also transcribed by its genome in the host cell. These viral RNAs have common 5ʹUTR and 3ʹUTR nucleotide sequences. The human miRNAs primarily target the 3ʹUTR of host mRNAs in the cytoplasm but may target the 3ʹUTR and 5ʹUTR of the infecting viral RNA genome and transcripts [23-26]. In the SARS-CoV-2 virus-infected human pancreas, the pancreas cell’s miRNAs may target the viral RNA UTRs instead of its native mRNAs due to the dominant presence of the viral RNA. This differential miRNA-targeting of the host genes may cause their upregulation in the cell.

We identified 21 human miRNAs (CoV-tar-miRNAs) potentially targeting the UTRs of the viral genome, among which 17 CoV-tar-miRNAs also target the Diabetes-associated genes, thus, regulating their expression before the infection. However, after infection, as the number of viral RNA copies increases in the cell, the SARS-CoV-2 genome may engage with the miRNAs regulating these native genes. This differential targeting of the miRNAs explains the upregulation of the Diabetes-associated genes after the viral infection (Fig. 2).

Fig. 2

Differential targeting of the host cell’s native miRNAs: a The cell’s native miRNAs regulate the genes by only targeting and suppressing the cell’s native mRNAs. b In a SARS-CoV-2-infected cell, the viral + ssRNA genome copies compete with the host cell’s native mRNAs in being targeted by the native miRNAs. The host miRNAs get apportioned in targeting the SARS-CoV-2 genome copies too. It reduces the contribution of miRNAs in regulating the host cell’s genes, thus, upregulating them

The SOCS3 (Suppressor of Cytokine Signaling 3) gene codes for a protein that helps regulate cytokine signal transduction [58]. Overexpression of SOCS3 has been observed in mice having Type 1 Diabetes. SOCS3-deficiency in pancreas beta cells is associated with increased resistance to apoptosis, thus, preventing Type 1 Diabetes [59]. We identified that the SOCS3 gene is the potential target of six 3ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-1236-3p, hsa-miR-3941, hsa-miR-4279, hsa-miR-4717-3p, hsa-miR-5088-5p and hsa-miR-6749-3p, and five 5ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-1303, hsa-miR-3123, hsa-miR-3925-5p, hsa-miR-4691-3p and hsa-miR-7851-3p. Among them, hsa-miR-3925-5p has been shown to have more favour towards the viral RNAs than the SOCS3 mRNAs. With the increase in the SARS-CoV-2 genome and transcript copies in the host pancreas cell, the viral RNAs may engage with the host miRNAs that were earlier targeting the SOCS3 mRNAs. Thus, SOCS3 upregulation after the SARS-CoV-2 infection may be due to the differential targeting of the host miRNAs.

The PSMB8 gene codes for Proteasome 20S Subunit Beta 8 and has been found to promote apoptosis [60, 61]. We see an enhanced expression of PSMB8 in the SARS-CoV-2-infected Pancreas tissue. The PSMB8 mRNA is the potential target of three 3ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-1236-3p, hsa-miR-4279 and hsa-miR-5088-5p, and three 5ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-298, hsa-miR-3123 and hsa-miR-3925-5p. These host miRNAs target the 3ʹUTR and 5ʹUTR of the SARS-CoV-2 genome. Among them, hsa-miR-298 and hsa-miR-3925-5p have been shown to have more chances of binding with the viral genome and transcripts than the PSMB8 mRNAs due to lower MFE with viral RNA. Thus, our study suggests that the upregulation of PSMB8 is caused due to differential targeting of the host miRNAs leading to apoptosis of the pancreas’ beta cells.

The CP gene encodes a secretory plasma protein called Ceruloplasmin, which is increased in the Diabetic condition [62]. We noticed that in the SARS-CoV-2-infected Pancreas tissue, the expression of the CP gene is enhanced. The CP mRNA is the potential target of two 3ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-466 and hsa-miR-4775, and four 5ʹUTR CoV-tar-miRNAs, i.e., hsa-miR-3123, hsa-miR-4691-3p, hsa-miR-4747-5p and hsa-miR-5196-5p. Among them, hsa-miR-4691-3p and hsa-miR-5196-5p have been shown to have lower MFE and more chances of binding to the viral genome and transcripts than the CP mRNAs. We suggest that as the number of viral RNA copies increases, the host microRNAs, i.e., hsa-miR-4747-5p and hsa-miR-5196-5p, may preferentially target the SARS-CoV-2 genome leading to upregulation of the CP gene. Our network analysis reveals that the highly connected CP gene has high closeness and betweenness centrality, indicating its strong influence on the network.

AGT gene encodes for the secretory pre-angiotensinogen or angiotensinogen precursor, which is an essential component of the renin-angiotensin system to maintain the blood pressure and fluid and electrolyte homeostasis in the body [63, 64]. The renin-angiotensin system’s expression in the pancreas has been shown to be enhanced in diabetic conditions [65]. The expression of AGT has been reported to be positively correlated with diabetes in rats [66]. Our study reveals that hsa-miR-1303, hsa-miR-298, has-miR-4645-3p, and hsa-miR-7851-3p are the 5ʹUTR CoV-tar-miRNAs, and has-miR-4775 is the 3ʹUTR CoV-tar-miRNA that target the AGT mRNAs. As the viral RNA copies increase in the cell, the differential targeting of these host miRNAs may lead to the upregulated expression of the AGT gene.

The CFB gene encodes for the secretory complement factor b, which links obesity to Diabetes. Its level is found to be increased during obesity and diabetes [67]. Obesity increases the risk of high COVID-19 severity. CFB is also linked to insulin resistance [68]. According to our study, the CFB mRNA is targeted by a 3ʹUTR CoV-tar-miRNA, i.e., hsa-miR-4775 and the differential targeting of this miRNA leads to its upregulation in the cell.

The MFE analysis revealed that hsa-miR-4691-3p and hsa-miR-5196-5p have lower MFE and higher chances of targeting the viral UTR region than the CP mRNA’s 3ʹUTR. Similarly, hsa-miR-3925-5p showed more chances of targeting the viral 5ʹUTR than the SOCS3 mRNA’s 3ʹUTR. hsa-miR-298 and hsa-miR-3925-5p showed more targeting probability with the viral 5ʹUTR than the PSMB8 mRNA’s 3ʹUTR. As per our knowledge, the expression of hsa-miR-466, hsa-miR-3123, hsa-miR-6749-3p, and hsa-miR-7851-3p in human pancreas cells is not known and must be investigated. The MFE analysis supports these four miRNAs to target the SARS-CoV-2 genome rather than the transcripts of Diabetes-associated DEGs. The MFE analysis indicates the possibility that AGT and CFB are not affected by the differential regulation by CoV-tar-miRNAs. However, the link of these genes to the associated CoV-tat-miRNAs suggests further experimental validation. The role of co-expression of CP, SOCS3 and PSMB8 on AGT and CFB must also be evaluated.

Among the Diabetes-associated CoV-tar-miRNAs, hsa-miR-466, hsa-miR-3123, hsa-miR-6749-3p and hsa-miR-7851-3p have not been confirmed to be expressed in the human pancreas. MIR-466 (the gene that encodes hsa-miR-466 in humans) has been reported to express in the rat’s pancreas [69]. However, the expression of hsa-miR-3123, hsa-miR-6749-3p and hsa-miR-7851-3p in the human pancreas is unknown and must be investigated in the future.

The literature search regarding the Diabetes-associated CoV-tar-miRNAs helped us understand their known relation with Diabetes and other pancreatic comorbidities. The Diabetes-associated CoV-tar-miRNAs, whether or not supported by MFE analysis, have been associated with Diabetes. Downregulation of hsa-miR-298 in the pancreas has been reported to contribute to mammalian pancreatic alpha cells’ resistance towards cytokine-induced apoptosis [70]. hsa-miR-298 has also been suggested as a diagnostic tool to predict Type 2 Diabetes [71, 72]. hsa-miR-3925-5p shows association with diabetic vascular complications induced by high glucose levels [73]. hsa-miR-4691-3p downregulation has been reported in the blood serum of individuals with Latent autoimmune diabetes in adults [74]. Its upregulation in adipocytes has been observed to be linked to obesity with increased insulin sensitivity following gastric bypass surgery [75, 76]. Leukocytes show downregulation of hsa-miR-5196-5p in gestational Diabetes mellitus [77].

hsa-miR-4717-3p, hsa-miR-6749-3p and hsa-miR-4645-3p show association with diabetic retinopathy [78, 79]. hsa-miR-4717-3p is also upregulated in periodontitis which is associated with diabetes [80]. hsa-miR-6749-3p is upregulated in pancreatic ductal adenocarcinoma and shows association with insulin resistance syndrome [81, 82]. hsa-mir-3123 is downregulated in insulin receptor haploinsufficient hepatic stellate cells [83]. hsa-miR-4775 is associated with diabetic nephropathy and its role in targeting apoptotic pathway’s core genes [84, 85]. hsa-miR-1236-3p is reported to be involved in apoptosis induction [86]. hsa-miR-1303 is upregulated in the blood serum of Type 2 Diabetes patients [87]. hsa-miR-7851-3p is downregulated in leukocytes during gestational Diabetes mellitus [88]. hsa-miR-466 is downregulated in Diabetic conditions and is associated with delayed wound healing in diabetes [72, 89–91]. hsa-miR-4279 and hsa-miR-3941 are associated with Diabetes. hsa-miR-4279 is expressed in urinary extracellular vesicles of Diabetic patients with macroalbumin, while hsa-miR-3941 is linked with islet immunity [92-94]. hsa-miR-5088-5p is upregulated in urine in the case of Diabetic Kidney disease [95]. hsa-miR-4747-5p is downregulated in human aortic vascular smooth muscle cells with high glucose-induced calcification/senescence [96]. It is also associated with insulin resistance of patients of post gastric bypass surgery obesity [97].

Overall, our study suggests that after SARS-CoV-2 infection, the transcripts of genes associated with Diabetes and cell death, i.e., CP, PSMB8, SOCS3, AGT, and CFB, compete with the high number of the viral transcripts in the host cell. The host miRNAs preferably target the viral genome copies and transcripts due to dominated presence in the host cell. Four miRNAs also show lower MFE with the viral RNAs than the host mRNAs. It indicates that their targeting of the viral genome and transcripts is similar to be what is expected and favored in nature. Thus, the Diabetes-associated genes get upregulated due to the differential miRNA targeting between the viral and host RNAs, which may also lead to the death of the host pancreas cell even if the miRNAs can tackle the infection by targeting the SARS-CoV-2 RNAs. The SARS-CoV-2 infection has been earlier reported to be associated with the induction of apoptosis in islet beta cells of the human pancreas [14]. Our study can also be applied to the insulin-producing pancreas beta cells, approving the low or no insulin secretion by them. This can worsen the condition of COVID-19 patients due to diabetic complications or the new onset of Diabetes, even if the pancreas cell’s miRNAs block the viral genome function. It can also lead to the hyperglycemic condition, reducing the individual’s immune response to handle the SARS-CoV-2 infection.

Our study suggests that artificial miRNAs can be designed and inserted into the infected cells for therapeutic purposes to bind with the viral genome, thus blocking both its function and the differential targeting of the host cell’s miRNAs. A research study with this concept has been done with four artificial miRNAs for targeting and blocking the Chikungunya viral genome [98]. The mechanism of differential miRNA targeting must be further validated in vitro.

Conclusion

This study highlights the effect of SARS-CoV-2-infection on the miRNA-regulation of Diabetes-associated genes in the human pancreas. The gene expression and disease enrichment analyses of the RNA-Seq data have revealed five biomarker genes, i.e., CP, SOCS3, AGT, PSMB8 and CFB, associated with Diabetes and get significantly upregulated in the pancreas following SARS-CoV-2-infection. In addition, we discovered four miRNAs, i.e., hsa-miR-298, hsa-miR-3925-5p, hsa-miR-4691-3p and hsa-miR-5196-5p, that may favorably target the SARS-CoV-2 RNAs over the Diabetes-associated upregulated gene transcripts. Thus, our study suggests that following the infection, the cell’s native miRNAs target the SARS-CoV-2 genome instead of the cell’s transcripts that were being targeted before the infection. This differential miRNA targeting causes the pancreas cell’s Diabetes-associated genes to upregulate, leading to diabetic complications or even new onset of Diabetes. Therefore, preventive, therapeutic methods are needed to block the viral genome from binding the host cell’s miRNAs and facilitate binding with the externally provided artificial miRNAs.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 5890 KB)