Competing endogenous RNA network profiling reveals novel host dependency factors required for MERS-CoV propagation

ABSTRACT Circular RNAs (circRNAs) are an integral component of the host competitive endogenous RNA (ceRNA) network. These noncoding RNAs are characterized by their unique splicing reactions to form covalently closed loop structures and play important RNA regulatory roles in cells. Recent studies showed that circRNA expressions were perturbed in viral infections and circRNAs might serve as potential antiviral targets. We investigated the host ceRNA network changes and biological relevance of circRNAs in human lung adenocarcinoma epithelial (Calu-3) cells infected with the highly pathogenic Middle East respiratory syndrome coronavirus (MERS-CoV). A total of ≥49337 putative circRNAs were predicted. Among the 7845 genes which generated putative circRNAs, 147 (1.9%) of them each generated ≥30 putative circRNAs and were involved in various biological, cellular, and metabolic processes, including viral infections. Differential expression (DE) analysis showed that the proportion of DE circRNAs significantly (P < 0.001) increased at 24 h-post infection. These DE circRNAs were clustered into 4 groups according to their time-course expression patterns and demonstrated inter-cluster and intra-cluster variations in the predicted functions of their host genes. Our comprehensive circRNA-miRNA-mRNA network identified 7 key DE circRNAs involved in various biological processes upon MERS-CoV infection. Specific siRNA knockdown of two selected DE circRNAs (circFNDC3B and circCNOT1) significantly reduced MERS-CoV load and their target mRNA expression which modulates various biological pathways, including the mitogen-activated protein kinase (MAPK) and ubiquitination pathways. These results provided novel insights into the ceRNA network perturbations, biological relevance of circRNAs, and potential host-targeting antiviral strategies for MERS-CoV infection.

Middle East respiratory syndrome coronavirus (MERS-CoV) is an emerging human-pathogenic coronavirus which has caused >2400 infections and >800 deaths worldwide since 2012 [17][18][19]. Despite a high mortality rate of >30%, the virus-host interaction and pathogenesis of this emerging infection remains incompletely understood [20]. Moreover, the host ceR-NAs profile changes and the functions of circRNAs in MERS-CoV and other human-pathogenic coronavirus infections have not been reported. In this study, we used MERS-CoV as a model to investigate the host ceRNA network changes and biological relevance of circRNAs in CoV infection. We showed that MERS-CoV significantly perturbed a high number of cir-cRNAs, miRNAs, and mRNAs which were involved in a wide range of biological processes. We also validated the effects of selected circRNAs in MERS-CoV replication and provided new insights into potential host-targeting antiviral strategies through the manipulation of circRNAs.

Materials and methods
Virus and cells MERS-CoV (strain HCoV-EMC/2012) was kindly provided by Ron Fouchier (Erasmus Medical Center, Rotterdam, the Netherlands) [21]. Clinical isolates of SARS-CoV (GZ50) and influenza A virus strain A/ HongKong/415742/2009(H1N1)pdm09 were obtained for validation studies from Department of Microbiology at The University of Hong Kong and prepared as previously described [22]. Calu-3 (human lung adenocarcinoma) cells were used to establish the MERS-CoV replication model for transcriptomic study as previously described according to Biosafety Level 3 practice [22,23]. The cells were maintained in Dulbecco's Modified Eagle Medium/F12 (DMEM/F12) supplemented with 10% heat-inactivated fetal bovine serum (FBS), 100 U/ml penicillin and 100 µg/ml streptomycin as previously described [24,25]. Primary human embryonic lung fibroblasts (HFL) were maintained in supplemented Minimum Essential Medium (MEM) as we described previously [26].

Sample preparation and total RNA isolation
Calu-3 cells (10 7 cells per biological replicate, three biological replicates), were mock-infected or infected with MERS-CoV at multiplicity of infection (MOI) of 4. The cell pellets were harvested at 24 h post-infection (hpi) in 1 ml of lysis buffer. Total RNA was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions, followed by DNase I digestion (Epicentre, Madison, WI, USA) for 15 min at 37°C. The integrity and quality of the extracted total RNA were evaluated using an Agilent 4200 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) with RNA integrity number (RIN) >7.0. The RNA quantity was measured using NanoDrop TM One (Thermo Fisher Scientific, Waltham, MA, USA).

Library construction and Illumina sequencing
A circRNA library, a small RNA library, and a mRNA library were constructed for the identification of cir-cRNA, miRNA, and mRNA, respectively. To prepare for circRNA sequencing, linear RNAs were removed using RNase R (RNR07250, Epicentre) (1 unit/ μg) for 20 min treatment at 37°C. Ribosomal RNA (rRNA) was depleted in the total RNA using a Ribo-Zero Gold Kit (Epicentre) following the manufacturer's instructions. After purification, the rRNA-depleted RNA products were fragmented using VAHTS Total RNA-seq (H/M/R) Library Prep Kit for Illumina (Vazyme Biotech Co., Ltd, Nanjing, Jiangsu, China). Three cDNA libraries were sequenced on Illumina Hiseq X-Ten platform (HaploX Biotechnology, Jiangxi, China) and 2×150 bp paired-end (PE150) reads were obtained using the HiSeq Control Software (HD3.5.0). Next, the sequencing reads were conducted for real-time sequencing image analysis and base-calling using Real-Time Analysis (v2.7.7). All Illumina sequencing raw and processed data were submitted to the GEO database (http://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE139516.

Identification of circRNAs, miRNAs, and mRNAs
The raw reads were subjected to quality assessment using fastp (0.19.5) [27]. Reads containing adapter, reads of N base over 5 bp, and low-quality reads were removed to obtain high quality clean reads. To identify circRNA, the clean reads were aligned with the human reference genome GRCh38 (http://genome.ucsc.edu/) using HISAT2 (https://bio.tools/hisat2). CIRI2 [28] and find_circ (version 1.2) [29] were used for the prediction of putative circRNAs. Overlapping circRNAs in CIRI2 and find_circ were selected for further analysis. The Burrows-Wheeler Alignment tool was used to identify miRNA [30]. Unique sequences containing 18-35 nucleotides were mapped to miRBase 22.0 by BLAST search to identify known and novel miRNAs. Due to the short lengths, the expression of circRNA and miRNA was normalized to transcripts per million (TPM), where TPM = (actual miRNA/circRNA counts of total clean read) × 10 6 . To identify mRNA, after alignment with the human reference genome GRCh38 with HISAT2, clean reads were quantified based on the number of reads spanning the back-splicing junction, and their fragments per kb for a million reads (FPKM) were calculated using HTSeq 0.10.0. RNAs with |log 2 (fold change)| ≥1 and adjust P value <0.05 were defined as differentially expressed (DE) by DESeq2 (version 1.18.1).

Weighted gene co-expression network analysis (WGCNA)
WGCNA was performed to identify mRNA-circRNA pairs with positive correlations. The WGCNA R package was downloaded from Bioconductor (https:// bioconductor.org/) and was applied to the dissimilarity matrix (1-TOM) to find clusters in each network. The DE mRNAs and DE circRNAs were clustered into different modules. To better understand the biological function of each module, we performed pathway analysis for all the DE circRNAs and DE mRNAs that belonged to the same module.
Integrated analysis of circRNA-miRNA-mRNA network circRNA-miRNA-mRNA interaction networks were constructed as described below. First, correlation analysis between DE circRNAs and DE miRNAs was performed and the correlation P value (cP) was calculated based on Pearson correlation coefficients. circRNA-miRNA pairs with strong negative correlations were defined as those having Pearson correlation coefficient (r) < −0.7 and cP <0.05. The circRNA-miRNA pairs with strong negative correlations were selected for circRNA-miRNA binding site prediction using miRanda (http://www. microrna.org/) v3.3a (-sc 140 -en −1.0 -scale 4 -out) [31]. Second, correlation analysis between DE miR-NAs and DE mRNAs was performed using the same method. The miRNA-mRNA pairs with strong negative correlations were again defined as those having r < −0.7 and cP value <0.05. The miRNA-mRNA pairs with strong correlations were selected for miRNA-mRNA binding site prediction using miRanda. Finally, correlation analysis between DE mRNAs and DE circRNAs was performed. The mRNA-circRNA pairs with strong positive correlations were defined as those having r >0.7 and cP <0.05. The final circRNA-miRNA-mRNA network graphs were constructed and visualized using Cytoscape v3.5.1 (http://www.cytoscape.org/).
RNase R resistance analysis of circRNAs and quantitative reverse transcription-polymerase chain reaction (qRT-PCR) The total RNA (1μg) extracted from the Calu-3 cells was treated with 1 unit of RNase R or nuclease-free water (mock control) and incubated for 20 min at 37°C. The digested RNA was purified using miRNeasy Mini Kit (Qiagen, Hilden, Germany). Divergent primers for the selected circRNAs were designed through CircInteractome (https://circinteractome.nia.nih.gov/) (Supplementary Table 1) [32]. Then, the treated RNAs were subjected to reverse transcription and quantitative polymerase chain reaction (qRT-PCR) with Transcriptor First Strand cDNA Synthesis Kit and Light-Cycler 480 master mix (Roche Holding AG, Basel, Switzerland) as previously described [26,33,34]. The housekeeping gene, GAPDH, was used as an internal control. The comparative CT (2 −ΔΔCT ) method was used to obtain the fold change of cir-cRNA expression levels.

siRNA knockdown and infection
Customized siRNAs were designed through CircInteractome and synthesized by Dharmacon (Lafayette, CO, USA). The detailed sequences of siRNAs were listed in Supplementary Table 1. Cells were transfected with 70nM siRNA using Lipofectamine RNAiMAX (Thermo Fisher Scientific) twice over two consecutive days as previously described [26,35]. At 24 h after the second transfection, the cells were challenged with MERS-CoV (MOI = 0.1). The inoculum was removed after 1 h at 37°C. The knockdown efficiency was assessed in parallel by qRT-PCR. At 24 hpi, the cells were harvested for further analysis. The cell viability of the siRNA-treated cells were evaluated with CellTi-terGlo® as described before [36].

Plasmids construction and transfection
Full lengths of hsa_circ_0067985 and hsa_-circ_0006275 were synthesized and subcloned into the pCD25-ciR vector and verified by sequencing in Geneseed Biotech Co.,Ltd (Guangzhou, China). The recombinant plasmids were transfected into Calu-3 cells using the Lipofectamine™3000 transfection reagent (Life Technologies) according to the manufacturer's protocol. The expression of each circRNA was measured at 36 h post-transfection to determine the transfection efficiency.

Statistical analyses
All data were analysed with GraphPad Prism (version 6.0, GraphPad, Inc) as we previously described [37]. The P-values were adjusted using Benjamini and Hochberg's approach to control the false discovery rate. Student's t-test was used to determine significant differences in gene expression changes in host cells after siRNA knockdown of selected circRNAs. Chisquare test was used to compare the number of DE cir-cRNAs, DE miRNAs, and DE mRNAs in MERS-CoV infection. One-way ANOVA was used to determine significant differences in the other parameters between different groups. P < 0.05 was considered statistically significant.

Landscape of circRNAs, miRNAs, and mRNAs identified in human lung epithelial cells
A schematic representation of the study design and computational analysis was shown in Figure 1. To optimize the conditions for RNA sequencing, we first determined the viral replication kinetics of MERS-CoV infection (MOI = 2 and 4) in Calu-3 cells by qRT-PCR. As shown in Figure 2A, there was ∼2-log increase and 4-log increase in viral load at 6 hpi (early phase) and 24 hpi (late phase). Next, to obtain a comprehensive landscape of the ceRNA network of the lung epithelial cells with either mock infection or MERS-CoV infection, we performed RNA sequencing on virus-infected Calu-3 cells (MOI = 4) to profile the expression patterns of circRNA, miRNA, and mRNA. To increase the confidence level of circRNA prediction, both CIRI2 and find_circ were used [28,29]. A total of 54605 and 99037 putative circRNAs were predicted by CIRI2 and find_circ, respectively ( Figure 2B). Among these putative circRNAs, 49337 were predicted by both CIRI2 and find_circ. These 49337 putative circRNAs had a significantly higher mean GC content than the miRNAs and mRNAs in our study ( Figure 2C). Among the 49337 putative circRNAs, 40961 (83.0%) were derived from coding regions ( Figure 2D). The majority (41690/49337, 84.5%) of the putative circRNAs containing a small number (1-5) of back splice reads ( Figure  2E). A minority (1986/49337, 4.0%) of them had head-to-tail junction with reads ≥ 20. The overall cir-cRNA expression profiles at these three time points were uniform and not significantly altered by MERS-CoV infection. Collectively, these basic characteristics of the putative circRNAs identified in our study were similar to those identified in other human and mammalian cells, including mouse P19 embryonic carcinoma, human neuroblastoma SH-SY5Y, and acute promyelocytic leukemia-derived NB4 [38][39][40].
The number of different circRNAs generated from individual host genes was highly variable ( Figure 2F). Among the 7845 genes which generated putative cir-cRNAs, 147 of them (1.9%) each generated ≥30 putative circRNAs (Supplementary Table 2). These 147 genes were shown to be involved in various biological, cellular, and metabolic processes in Gene Ontology (GO) enrichment analysis (Supplementary Figure 1). Moreover, many of these genes were involved in viral infection. For example, protein tyrosine phosphatase receptor kappa (PTPRK) is downregulated by Epstein-Barr virus in Hodgkin lymphoma (HL), leading to suppressed transforming growth factor-beta (TGF-β) signalling [41]. Membrane-associated guanylate kinase with inverted organization-1 (MAGI-1) is a PDZ protein that interacts with viral proteins containing the PDZ-ligand binding motif, such as influenza A virus NS1, human papillomavirus 16 E6, and human alphapapillomavirus E6 proteins [41][42][43]. Notably, the different putative circRNAs generated from the same host gene might demonstrate diversified expressions. For example, the expressions of the different putative circRNAs generated from DNA polymerase epsilon subunit 2 (POLE2) and Ubiquitin Specific Peptidase 48 (USP48) varied by >2 folds at 24 hpi ( Figure 2G). These diversifications in the expression of the individual circRNAs might be related to their different biogenesis mechanisms and biological functions. circRNAs may be positively correlated, negatively correlated, or non-correlated with their corresponding host genes [4,29,44,45]. In MERS-CoV infection, we found that 922/1815 (50.8%) and 190/1815 (10.5%) of the DE circRNAs were significantly positively and negatively correlated with their host genes, respectively ( Figure 3C). To investigate the potential functions of these DE circRNAs, we performed pathway and functional enrichment analyses on their host genes using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and GO analysis, respectively. As shown in Figures 3D and 3E, a wide range of pathways and processes were perturbed. The 10 most upregulated DE circRNAs at 24 hpi and the functions of their host genes were illustrated in Figure 3F. These virus-perturbed functions included ATP metabolism, protein transport, protein phosphorylation and polyubiquitination, mitosis, regulation of translational initiation, and autophagy.

MERS-CoV infection induced global changes in
The DE circRNAs in different clusters exhibited inter-cluster and intra-cluster variations in host pathway perturbations Next, we conducted K-means clustering analysis to determine the expression kinetics of the DE circRNAs. As shown in Figure 4A, the point of inflexion fell between 3 and 5 in the Elbow method which prompted us to cluster the DE circRNAs into 4 groups (Clu1, Clu2, Clu3, and Clu4). These DE circRNAs demonstrated different time-course expression patterns ( Figure 4B). The number of DE circRNAs in each cluster was 1046, 239, 347, and 183, respectively, for Clu1 to Clu4. The mean expression intensities of the DE cir-cRNAs in Clu1, Clu2, and Clu3 in mock-infected cells were similar to those in MERS-CoV-infected cells at 6 hpi ( Figures 4B). In contrast, the mean expression intensity of the DE circRNAs in Clu4 in MERS-CoVinfected cells at 6 hpi was lower than that of mockinfected cells. The DE circRNAs in all 4 clusters in MERS-CoV-infected cells at 24 hpi were higher than those in mock-infected cells and MERS-CoV-infected cells at 6 hpi. Using KEGG pathway analysis, we showed that the predicted functions of the host genes of the DE circRNAs in the 4 clusters demonstrated inter-cluster and intra-cluster variations ( Figure 4C). Similar to the KEGG analysis of the host gens of the total DE circRNAs ( Figure 3D), ubiquitin-mediated proteolysis was the most perturbed pathway in both Clu3 and Clu4 ( Figure 4C). This observation implicated that key DE circRNAs that were involved in the regulation of MERS-CoV infection might be enriched in these clusters.
Ago2 protein is a core component of RNA-induced silencing complex that binds miRNAs to target mRNAs [47]. Explorations on the RBPs that could bind to these circRNAs through CircInteractome showed that both circFNDC3B and circCNOT1 contained multiple binding sites for Ago2 protein and supported our postulation that these circRNAs could act as miRNA sponges ( Figure 7A) [32]. We designed divergent qRT-PCR primers and performed RNase R resistance experiments to examine the expression of these circRNAs with or without MERS-CoV infection (Supplementary Table 1). As shown in Figure 7B, no significant reduction in the relative expression of the two circRNAs was detected after RNase R treatment, while the expression of the mRNAs were greatly affected. To examine the biological relevance of these circRNAs during MERS-CoV infection, siRNAs targeting the back-splicing site of circFNDC3B and cir-cCNOT1 were synthesized and transfected into Calu-3 cells to assess whether depletion of these two candidate circRNAs could limit virus replication (Supplementary Table 1). Our results confirmed that our siRNAs could efficiently knock down the circRNA expression ( Figure 7C) but not the linear form (Supplementary Figure 2), and most of the tested siRNAs significantly reduced MERS-CoV (but not SARS-CoV or influenza A/H1N1 viruses) replication in both cell lysate and culture supernatant ( Figure 7D and Supplementary Figure 3). Importantly, the siRNAs did not affect cell viability (Supplementary Figure 4). Upon successful knock down (Supplementary Figure  5A), the viral load reduction by si-circFNDC3 and si-circCNOT1 were similarly observed in HFL cells Figure 6. Potential viral pathogenic circRNAs in the ceRNA co-regulatory network. circRNAs, miRNAs, and mRNAs potentially involved in MERS-CoV pathogenesis were represented by triangle, rectangle, and circle nodes, respectively. The border thickness and filling colour of each node were mapped according to the expression and adjust P value of each RNA at 24 h post MERS-CoV infection. The size of mRNAs was proportional to their correlation extent with selected circRNAs. The stronger that correlation, the larger the node. Top 10 overrepresented GO terms were adopted to colour the border of mRNAs identified in the interactome, and the mRNAs strongly correlated with miRNAs (r < −0.85, cP < 0.05) and circRNAs (r > 0.85, cP < 0.05) simultaneously were labelled with name. Edge thickness was proportionally correlated with the predicted interaction between each circRNA-miRNA pair and miRNA-mRNA pair as defined by miRanda. Among the 6 circRNAs which were significantly upregulated, circ_0006275 and circ_0067985 were selected for further validation. The expression levels of circ_0032503, and the target miRNA of circ_0001680 and circ_0001524 were low, and were therefore not suitable for siRNA knockdown experiment. circ_0029617 only interacted with 1 miRNA and was therefore not selected for further validation. Figure 5B). Moreover, over-expression of circFNDC3B and circCNOT1 in Calu-3 cells enhanced virus replication (Supplementary Figure 6). To validate our hypothesis that circFNDC3B and cir-cCNOT1 function by sequestering their target miRNAs to regulate mRNA expression, we evaluated the expression of their representative target mRNAs after knockdown with si-circFNDC3B-2 or si-circCNOT1-1 which were selected based on their higher inhibitory activity of MERS-CoV replication over si-circFNDC3B-1 or si-circCNOT1-2, respectively. As shown in Figure 7E, circFNDC3B depletion significantly reduced the expression of its representative target mRNAs, including MAP3K9, MYO15B, and SPOCK1. Similarly, circCNOT1 depletion reduced the expression of its target mRNAs, including MAP3K9, MEF2C, USP15, and ZBTB11. Taken together, these results supported our hypothesis that the circRNAs identified in our circRNA-miRNA-mRNA network significantly impacted MERS-CoV replication.

Discussion
As an integral component of the ceRNA network, cir-cRNAs harbouring MREs can regulate gene expression by functioning as miRNA sponges to release the inhibitory effects of miRNAs on target genes [6,7]. The sponging effect of circRNAs has been shown to be more efficient than that of the linear miRNA and lncRNA transcripts [29,48]. It has previously been shown that various DNA and RNA viruses (herpesviruses, polyomaviruses, retroviruses, and hepaciviruses) could modulate viral replication via miRNA- (C) qRT-PCR examining the knockdown effect of siRNA candidates. Linear RNA of GAPDH was used as internal reference for normalization. *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001, one-way ANOVA. (D) Depletion of circFNDC3B and circCNOT1 suppressed MERS-CoV replication in cell lysate and supernatant. Scramble siRNA was served as a negative control. *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001, one-way ANOVA. (E) CircFNDC3B and circCNOT1 knockdown decreased the expression of representative targeting genes. Student's t-test was adopted to calculate the significance of gene expression with or without MERS-CoV infection. *P-value < 0.05; **P-value < 0.01; ***P-value < 0.001; ****P-value < 0.0001. mediated gene regulation [49,50]. In contrast, the role of circRNA in viral replication is less well understood. For coronaviruses, only two circRNA profile analyses on the animal-pathogenic PEDV and TGEV have been reported [13,15]. In this study, we used the highly human-pathogenic MERS-CoV as a model to demonstrate the interactions of circRNAs with other major components of the host cell ceRNA network and validated the effects of these circRNAs on coronavirus replication.
In our comprehensive profiling of the circRNAs, miRNAs, and mRNAs in human lung epithelial cells with or without MERS-CoV infection, a number of important observations were made. First, a large number of circRNAs (49337) were predicted, including 16285 (33.0%) previously annotated circRNAs and 33052 (67.0%) novel putative circRNAs. Most of these circRNAs were derived from coding regions (83.0%). This is in line with the findings from PEDVinfected porcine intestinal epithelial (IPEC-J2) cells [13], and may represent a conserved character of coronavirus-infected cells. Second, most (84.5%) of these putative circRNAs had low expression levels as evidenced by their small number (1-5) of supported reads. This finding corroborated with those identified in PEDV-infected porcine intestinal epithelial cells, avian leukosis virus subgroup-J-infected chicken hepatic cells, and simian vacuolating virus 40-infected Vero cells [8,10,13]. Third, the mean GC content percentage of the circRNAs was significantly higher than those of the miRNAs and mRNAs. As the GC content is positively correlated with the stability of RNA transcripts, our finding corroborated with the consensus that cir-cRNAs are usually highly stable [29,38]. Finally, our data showed that the number of different circRNAs generated from individual host genes was highly variable, suggesting a differential potential of gene regulation by diverse host genes through circRNA generation.
Our comprehensive circRNA-miRNA-mRNA network showed that MERS-CoV induced significant changes in the expression of many host cell circRNAs, miRNAs, and mRNAs. Compared with the circRNA-miRNA-mRNA network in Madin-Darby Canine Kidney (MDCK) cells infected with influenza A/H3N2, our network in MERS-CoV-infected Calu-3 cells demonstrated a higher number of DE circRNAs (7 vs 3), miR-NAs (19 vs 1), and mRNAs (547 vs 9). This might imply that MERS-CoV infection induced a more profound and global change in the host ceRNA network compared to the less virulent influenza A/H3N2, although the differences in experimental set up should also be considered [14]. The DE circRNAs in our network were associated with a wide range of biological, cellular, and molecular processes. Interestingly, using both KEGG pathway and GO functional analyses, we showed that ubiquitin-mediated proteolysis was significantly perturbed in MERS-CoV infection. The papain-like protease of MERS-CoV exhibits deubiquitinating activity and is involved in proteolysis of the viral polyprotein during virus replication [51,52]. Inhibitors of MERS-CoV papain-like protease such as 6-mercaptopurine and 6-thioguanine exhibit antiviral activity in vitro [18]. Modulation of the circRNAs associated with ubiquitin-mediated proteolysis identified in our study may provide a new antiviral strategy for MERS-CoV infection.
To validate the biological relevance of the DE cir-cRNAs identified in our network, we selected two DE circRNAs and investigated their effects on MERS-CoV replication and the expression of their target genes in human lung epithelial cells with or without siRNA knockdown. Our results showed that specific knockdown of circFNDC3B and circCNOT1 significantly reduced MERS-CoV viral load in both Calu-3 and HFL cells, which was potentially associated with the downregulation of circFNDC3B-and circCNOT1regulated target genes. For example, MAP3K9 is an upstream modulator of the mitogen-activated protein kinase (MAPK) pathways which influences many aspects of cell proliferation, migration, and apoptosis [53]. The extracellular signal-regulated kinase (ERK)/ MAPK signalling response is specifically modulated in MERS-CoV infection [54]. In this regard, we showed that siRNA knockdown of either circFND3B or cir-cCNOT1 resulted in significantly reduced expression of MAP3K9 and provided novel insights to modulate the ERK/MAPK pathway as a host-targeting antiviral strategy for MERS-CoV infection.
In addition to MAP3K9, siRNA knockdown of circFNDC3B or circCNOT1 similarly resulted in significantly reduced expression of other target genes. Ubiquitin-specific protease 15 (USP15) promotes RIG-I-mediated antiviral signalling by interacting with ubiquitin E3 ligase tripartite motif protein 25 (TRIM25) and plays important roles in the replication of various DNA and RNA viruses, including human papillomaviruse (HPV), human immunodeficiency virus (HIV), and hepatitis C virus (HCV) [55][56][57][58]. Myocyte enhancer factor 2C (MEF2C) is a transcription factor that is associated with the super-enhancer activity and cancer progression of Epstein-Barr virus infection [59]. The functions of SPARC/Osteonectin, Cwcv and Kazal-like Domains Proteoglycan 1 (SPOCK1), Zinc Finger and BTB Domain Containing 11 (ZBTB11), and myosin XVB (MYO15B) in viral infection are unclear.
In summary, our study provided novel insights into the ceRNA network perturbations and biological relevance of circRNAs in MERS-CoV infection. Knockdown of specific DE circRNAs in MERS-CoV infection resulted in significantly reduced viral load and may pave new ways for host-targeting antiviral strategies for this highly virulent emerging virus.