Effect of SARS-CoV-2 infection on host competing endogenous RNA and miRNA network

Data

Single-cell RNA-Seq and bulk miRNA-Seq data of Severe acute respiratory syndrome-related coronavirus (SARS-CoV-1) and Severe acute respiratory syndrome coronavirus (SARS-CoV-2) infected human cell lines were downloaded from Single Cell Expression Atlas (accession number PRJNA625518) (Wyler et al., 2021; Papatheodorou et al., 2019) and Gene Expression Omnibus (GEO accession no GSE148729) (Barrett et al., 2012), respectively. miRNA-Seq reading data had been normalized to reads per million. To construct miRNA:target interaction network, strongly supported miRNA and gene target pairs were downloaded from miRTarBase (version 8) (Huang et al., 2019) which houses experimental miRNA:target pair data. Human miRNAs targeting SARS-CoV-2 genome were extracted from four previous reports which utilized various prediction algorithms (Khan et al., 2020; Chow & Salmena, 2020; Demirci & Adan, 2020; Fulzele et al., 2020). The miRNAs common in at least three studies were considered as targeting the SARS-CoV-2 genome. SARS-CoV-1 targeting miRNAs were extracted from the study of Khan et al. (2020).

As viral RNA read counts, we considered the number of mapped reads (per million reads) of viral RNAs (i.e., SARS-CoV-1/2) at 12 h post-infection in Calu3 cells (Wyler et al., 2021) to simulate its competitive behavior with human gene transcripts. Viral RNA read counts were normalized to viral genome length in kb (i.e., RPKM). The single-cell RNA-Seq data and bulk miRNA-Seq data from SARS-CoV-1/2 infected Calu3 cells were combined according to miRNA:target pair dataset from miRTarBase (Huang et al., 2019). The perturbation simulations were carried out with and without viral RNAs.

ceRNA network simulations in single-cell datasets

The ceRNA network simulations were carried out using the ceRNAnetsim R package we had developed previously (Yuka & Yilmaz, 2021; Yuka & Yilmaz, 2020). Briefly, each node in the network was used as a trigger for perturbation after which perturbation ability in the whole network was calculated. Node importance was evaluated by calculating the average change in expression for all nodes in miRNA:target network when that node was perturbed. Simulation with 2,800 cells for each group (i.e., SARS-CoV-1, SARS-CoV-2, and Mock samples) that were sampled, was performed with 20 iterations after 3 fold up-regulation of each node in parallel (Tange, 2011). SARS-CoV-1/2 dataset simulations were performed in miRNA:target network with/without viral RNAs. Viral RNA was incorporated into the ceRNA network through miRNAs targeting the viral RNA.

Detecting miRNA:target interactions disturbed in presence of viral RNAs

We performed perturbation simulation using expression data for each cell in the ceRNA network in the absence and presence of viral RNA and calculated the perturbation ability of each node. Threshold value to be considered perturbing node was calculated by finite mixture model in results of all cell types (mock, SARS-CoV-1 infected, SARS-CoV-2 infected) (Trang et al., 2015). By comparing the perturbation ability of nodes in the presence and absence of SARS-CoV-1/2 RNA, we revealed nodes with altered perturbation ability in presence of SARS-CoV-1 or SARS-CoV-2.

In order to find nodes that are affected by genes with altered perturbation ability, we performed simulations for each gene with cells in which the gene lost its perturbation. Simulations of genes that lost perturbation ability in SARS-CoV-2 infection were carried out under the same conditions (20 iterations after 3 fold up-regulation) to find the genes affected by perturbation of these nodes. Functional annotation analysis was performed with DAVID and TAM 2.0 functional annotation tools were for genes and miRNAs, respectively (Jiao et al., 2012; Li et al., 2018).

Single-cell expression data overlaid on ceRNA network originated from miRTarBase

Single-cell RNA (sc-RNA) sequencing data of Calu3 cells at 12 h after infection contained 25,757 gene transcripts for 3,779 SARS-CoV-1 infected cells and 14,283 SARS-CoV-2 infected cells. In addition, sc-RNA sequencing data contained 10,501 mock cells at 12 h. Bulk reading counts for 1,541 miRNAs were available at 12 h from small RNA-Seq data. After integrating all the data, a ceRNA network comprised of 9,009 miRNA:target interactions among 739 miRNAs and 3,032 targets were constructed (Table 1). For each single cell, the expression data was overlaid on the ceRNA network prior to simulations.

For ceRNA networks containing SARS-CoV RNA, existing networks were extended with miRNAs potentially target viral RNAs (Khan et al., 2020; Chow & Salmena, 2020; Demirci & Adan, 2020; Fulzele et al., 2020). 129 miRNAs targeted SARS-CoV-1 RNA, 39 of which have already been in the ceRNA network, so the network was extended by 90 miRNAs. For SARS-CoV-2, 51 miRNAs targeted SARS-CoV-1 RNA, 14 of which have already been in the ceRNA network, so the network was extended by 37 miRNAs (Fig. S1). The contents of extended ceRNA networks are summarized in Table 1.

Since mock cells lacked miRNA reads at 12 h, we performed differential expression analysis using DEseq2 package with miRNA reads (TPM) at 4 and 24 h (Love, Huber & Anders, 2014). We found that 13 (miR-135b-5p, miR-141-3p, miR-186-5p, miR-18a-5p, miR-18b-5p, miR-19a-3p, miR-20a-5p, miR-29a-3p, miR-29c-3p, miR-30e−5p, miR-424-5p, miR-539-3p, miR-582-5p) of 1541 miRNAs were differentially expressed between 4 and 24 h in Mock cells (Fig. S2). Since the number of differentially expressed miRNAs was negligible, we considered 24-hour mock count data as comparable with 12-hour infection data.

Perturbation ability of miRNAs and genes in calu3 cells changes upon infection

Simulations in 2,800 of each type of cells (i.e., mock, CoV-1- and CoV-2-infected) revealed that particular miRNA and gene nodes in ceRNA networks had altered perturbation abilities in the absence and presence of viral RNA. Since the simulations were performed in many cells, the finite mixture model was used to determine the minimum number of cells required for a gene or miRNA node to be considered having perturbation ability. Since the perturbation simulations were performed in many cells, a cut-off, the minimum number of cells, was calculated above which the node is considered significantly perturbing. Perturbing miRNA and target numbers indicate the number of miRNAs and target genes which are above the corresponding cut-off value (Table 2).

We performed ceRNA network simulations of Mock, and infected cells with and without viral RNAs (CoV-1 and CoV-2). When we compared perturbation simulation results of three ceRNA networks (Mock, Infected without viral RNA, Infected with viral RNA) for each SARS-CoV, pairwise Pearson correlations of Mock vs w/o SARS-CoV-1 RNA and Mock vs w/o SARS-CoV-2 were calculated as 0.96 and 0.91, respectively (Figs. S3–S4). Thus, analysis of changed perturbation ability was performed on results from networks with and without viral RNA only.

As a result of perturbation simulations in the ceRNA network of SARS-CoV-1 infected cells, it was observed that SARS-CoV-1 RNA was perturbing in almost all cells (2779 of 2800), while 6 genes (HMBOX1, MYCN, IGFBP3, DIRAS3, IKZF1, ZNF763) lost their perturbation ability in presence of SARS-CoV-1 RNA in ceRNA network (Supplementary Data). Interestingly, 7 miRNAs (miR-1229-3p, miR-1307-3p, miR-182-5p, miR-22-5p, miR-30c-1-3p, miR-425-3p, miR-550a-5p) gained perturbation ability after SARS-CoV-1 integration (Table 3).

On the other hand, trends of changes in perturbation ability in SARS-CoV-2 infected cells were drastically different. It was observed that 195 genes lost their perturbation ability in the SARS-CoV-2 RNA integrated ceRNA network simulations (Supplementary Data). In contrast to SARS-CoV-1 results, 29 miRNAs gained perturbation ability in presence of SARS-CoV-2 RNA (Table 3). Please note that no miRNA with lost perturbation ability was detected in simulations in presence of both viral RNAs. Only one gene and 2 miRNAs were found to be common in the list of nodes that had altered perturbation ability in presence of each SARS-CoV RNAs. That gene and two miRNAs were listed in Tables S1–S2.

Functional analysis of genes losing perturbation ability only in SARS-CoV-2 infection

A total of 195 genes that lost perturbation ability specifically in presence of SARS-CoV-2 RNA were subject to functional analysis. The genes which lost perturbation ability showed enrichment in several important KEGG pathways such as p53 signaling pathway (hsa04115, 8.48 fold enrichment with 0.000002 p-value), Hepatitis B (hsa05161, 4.31 fold enrichment with 0.0002 p-value), and Cytokine-cytokine receptor interaction (hsa04060, 3.27 fold enrichment with 0.0003 p-value). We observed that they enriched in many molecular functions as transmembrane receptor protein serine/threonine kinase activity (GO:0004675, 41.91 fold enrichment with 0.0001 p-value), protein kinase activity (GO:0004672, 3.415026 fold enrichment with 0.0004 p-value), and activin binding (GO:0048185, 39.294693 fold enrichment with 0.000005 p-value). Also, in important biological processes such as negative regulation of ERK1 and ERK2 cascade (GO:0070373, 12.38 fold enrichment with 0.000003 p-value), positive regulation of protein kinase B signaling (GO:0051897, 8.55 fold enrichment with 0.00004 p-value), and response to drug (GO:0042493, 4.73 fold enrichment with 0.000001 p-value) significant enrichment were determined. All data is available in the Supplementary Data.

In order to decipher the potential downstream effects of SARS-CoV-2 RNA in the ceRNA network, we identified genes that are affected through the ceRNA network by 195 genes (Fig. 1A). First, we collected barcodes of cells where perturbation was lost for each of the 195 genes in the presence of SARS-CoV-2 RNA. Then we performed simulations triggered from each gene and found the most affected genes from perturbation simulations. For example, in simulations using 102 cells for the ASVR1B gene, CCND1 gene was affected in 91 cells (asterisk, ratio 0.89). In the simulations performed with the CD47 gene in 103 cells, we found that CCND1 gene was among the most affected genes in 56 cells (triangle, ratio 0.54) (Fig. 1A). Ultimately, 18 genes (TNFRSF10B, STAT3, PKM, MYH9, MYC, MAP4K4, IL6, HMGA1, HMGA2, E2H2, EWSR1, DICER1, CDK6, CCND1, CASP8, APP, AGO1, ACTG1) appear to be affected by almost all 195 genes in absence of SARS-CoV-2 RNA and were no longer perturbed in presence of SARS-CoV-2 RNA.

Functional annotation analysis of the affected genes via DAVID (Jiao et al., 2012) showed enrichment in several important KEGG pathways (Fig. 1B). For instance, Hepatitis B (hsa05161, 18.98 fold enrichment with 0.000007 p-value), measles (hsa05162, 17.24 fold enrichment with 0.0001 p-value), microRNAs in cancer (hsa05206, 9.62 fold enrichment with 0.0001 p-value), and viral carcinogenesis (hsa05203, 11.18 fold enrichment with 0.0006 p-value) were among pathways with highest enrichment. Interestingly, these genes were significantly enriched in biological processes such as response to estradiol (GO:0032355, 43.42 fold enrichment 0.00008 p-value) and oncogene-induced cell senescence (GO:0090402, 987.76 fold enrichment with 0.001 p-value). All of the genes were also enriched significantly in molecular functions such as protein binding (GO:0005515, 1.92 fold enrichment with 0.00002 p-value) and transcription factor binding(GO:0008134, 17.48 fold enrichment with 0.0001 p-value). Full functional annotation data of affected genes is available in the Supplementary Data.

Our analysis revealed that specifically in presence of SARS-CoV-2 RNA, several genes lose their perturbation ability in the ceRNA network, rendering particular genes uncontrollable by the ceRNA network which were shown to be associated with diseases or cellular functions in literature. We found three genes (CASP8, ACTG1, and CCND1) that were significantly associated with viral myocarditis. Although several reports highlighted the issue of SARS-CoV-2-induced viral myocarditis, evidence for the involvement of these genes in the onset of viral myocarditis is lacking (Siripanthong et al., 2020). Although its molecular mechanisms have not been explained yet, studies on COVID-19 patient samples had reported the loss of thyroid function (CCND1, MYC, ACTG1 genes in Thyroid hormone signaling pathway with 11.96 fold enrichment and 0.02 p-value, was observed) and related loss of immune system regulation (Lui et al., 2020; Khoo et al., 2021). SARS-CoV-2-induced viral carcinogenesis has not been addressed thoroughly yet but in a recent study, the effect of the virus on tumor development had been investigated using SARS-CoV-2-host networks and the potential role of virus for tumorigenesis had been suggested (Souchelnytskyi, Nera & Souchelnytskyi, 2021). Interestingly, we found that 4 genes (CASP8, CCND1, MYC, and STAT3) were enriched in the biological process of response to estradiol. In the literature, it had been reported that estradiol may be a therapeutic option in SARS-CoV-2 viral infection due to the immune regulatory function (Mauvais-Jarvis, Klein & Levin, 2020). For detailed functional annotation and enrichment results please refer to the Supplementary Data. Although there are many recent studies on the systemic effects or symptoms of SARS-CoV-2 infection, number of studies that dealt with SARS-CoV-2-induced regulatory mechanisms on a large scale network is quite limited.

Functional annotations of miRNAs that gained perturbation ability in the presence of SARS-CoV-2 RNA

We observed that 29 miRNAs gained perturbation ability as a result of the integration of SARS-CoV-2 RNA to the host ceRNA network. By using TAM2.0 database (Li et al., 2018), we observed that miRNAs were significantly enriched in a wide variety of diseases such as metabolic disorders and carcinoma. They were also enriched in many cellular functions such as inflammation and immune response (Fig. 2).

Strikingly, we observed that several functions and diseases associated with miRNAs overlapped with symptoms of COVID-19 cases. For example, miR-23b, miR-30c-1-3p, miR-424-3p, miR-22, miR-146a, let-7e, miR-539, and miR-29a were found to be associated with heart failure (4.08 fold enrichment with 0.0004 p-value) which accompanies SARS-CoV-2 infections (Tomasoni et al., 2020). miR-22 and miR-146a were significantly enriched in emphysema, shortness of breath, (61.28 fold enrichment with 0.0002 p-value) which was reported as a strong complication of SARS-CoV-2 infection (Lacroix et al., 2020). Previous reports had shown that miR-22 controls the activation of immune system cells through the activation of histone deacetylase HDAC4, while miR-146a promotes the generation of abnormal inflammatory responses through a fibroblast-mediated way (Lu et al., 2015; Sato et al., 2010). These miRNAs may contribute to the development of chronic inflammatory disease through different mechanisms. It was reported that SARS-CoV-2 infection had the potential to cause Diabetes Mellitus, Type 2 through damage to β-cells (Hayden, 2020) and several miRNAs (miR-766, miR-652, miR-29a, miR-409, miR-23b, miR-15a, miR-22, miR-146a, let-7e) were significantly enriched in this disease. Our results pointed out direct or indirect associations between COVID-19 related diseases and the host ceRNA network.

Although there is not full knowledge about the molecular mechanisms regarding destructive effects of COVID-19, miRNA interactions of viral infection may have significant contributions. Considering the miRNA:ceRNA interactions in SARS-CoV-2 infections may facilitate the understanding of the systemic dysfunctions caused by viral RNA.