Profiling and functional analysis of differentially expressed circular RNAs identified in foot-and-mouth disease virus infected PK-15 cells

Circular RNAs (circRNAs) are a new type of endogenous noncoding RNA that exhibit a variety of biological functions. However, it is not clear whether they are involved in foot-and-mouth disease virus (FMDV) infection and host response. In this study, we established circRNA expression profiles in FMDV-infected PK-15 cells using RNA-seq (RNA-sequencing) technology analysis. The biological function of the differentially expressed circRNAs was determined by protein interaction network, Gene Ontology (GO), and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway enrichment. We found 1100 differentially expressed circRNAs (675 downregulated and 425 upregulated) which were involved in various biological processes such as protein ubiquitination modification, cell cycle regulation, RNA transport, and autophagy. We also found that circRNAs identified after FMDV infection may be involved in the host cell immune response. RNA-Seq results were validated by circRNAs qRT-PCR. In this study, we analyzed for the first time circRNAs expression profile and the biological function of these genes after FMDV infection of host cells. The results provide new insights into the interactions between FMDV and host cells.


Introduction
Circular RNAs (circRNAs) belong to a class of noncoding RNAs that are widespread in the cytoplasm of eukaryotic cells and are structurally and functionally different from linear RNA molecules [1]. They are covalently closed-loop RNA molecules that are formed by back-splicing of the 5′ and 3′ ends of the primary transcript. They have a strong structural stability, tissue, and spatiotemporal specificity [2,3]. To date, three types of circRNA molecules have been reported including cir-cRNAs generated by reverse exon splicing, circRNAs that form by intronic lasso, and circRNAs consisting of both introns and exons [4][5][6][7]. Previously, circRNAs were considered by-products of abnormal splicing during transcription; however, with the rapid development of high-throughput sequencing technologies and bioinformatics, there is growing evidence that circRNAs are involved in regulating a variety of important physiological functions [8]. Most circRNAs are composed of exons and are located in the cytoplasm, indicating that they function as protein regulators in the translation and modification of proteins [7,[9][10][11][12]. Recently, it has been shown that circRNAs act as sponges for microRNAs (miRNAs), which act as competing endogenous RNAs (ceRNAs) to regulate post-transcriptional gene expression events [13][14][15]. Other studies have found that cir-cRNAs play important regulatory roles in pathological processes such as neurological diseases [16] and cancer [17,18]. In addition, because of the low molecular weight of circRNAs, they are transported by extracellular vesicles, such as nanoparticles and exosomes. They are now widely studied and considered as molecular markers, therapeutic targets, and drug carriers in a variety of diseases [19].
As novel regulatory molecules, circRNAs mediate the regulation of viral infections and the cellular immune response, which provide a new perspective for understanding virus-cell interactions. Recently, several studies have shown that circRNAs are also involved in the regulation of virus-host cell interactions as well as in the antiviral cell immune response [5,20]. For example, during H1N1 influenza A virus (IAV) infection, circGATAD2A overexpression in the host promotes replication of the H1N1 IAV virus [21]. Hepatitis C virus infection induces host circRNAs to exert nonsense-mediated decay and inhibit viral replication [22]. Viruses can also use circR-NAs to interfere with the host antiviral immune response and help to escape immune surveillance and antiviral immunity. Interestingly, when viral invasion occurs, some circRNAs inhibit immune cell activation; however, when circRNAs are degraded by RNase L, they subsequently regulate autoimmune disease or viral infection clearance by activating PKR activation and downstream cascade responses [23]. Although studies have demonstrated that circRNAs are involved in the regulation of host cells after viral infection, there is a lot to unravel to understand the role that circRNAs play in virus-host interactions and viral pathogenesis.
Foot-and-mouth disease (FMD) is a highly contagious viral disease that occurs in livestock worldwide. FMD is caused by the foot-and-mouth disease virus (FMDV), which mainly infects cloven-hooved animals [24]. Pathological blisters appear in the oral mucosa, extremities, and breasts as the main clinical symptoms, and animals may die from severe infections. This causes significant losses to the animal husbandry industry and the economy [25]. FMDV has a genome of approximately 8.5 kb and is a single-stranded positive-sense RNA virus of the genus Aphthovirus, within the family Picornaviridae. It contains an open reading frame (ORF) that encodes four structural proteins (VP1, VP2, VP3, and VP4) and 10 nonstructural proteins (L pro , 2A, 2B, 2C, 3A, 3B1, 3B2, 3B3, 3C, and 3D) [26][27][28]. FMDV is divided into seven serotypes (A, O, C, SAT 1, SAT 2, SAT 3, and Asia 1) according to geographical distribution and each serotype has a wide range of antigenic characteristics which do not elicit effective cross-protection. The high incidence and extensiveness worldwide make the prevention and control of FMD a great challenge [29]. Several studies have shown that FMDV is able to escape the host innate immune response through multiple pathways. For example, FMDV antagonizes host antiviral interferon production by inducing PERK and AKT-MTOR signaling to regulate autophagy [30,31]. FMDV L pro is able to target ISG15 by specifically cleaving the peptide bonds before the C-terminal Gly-Gly motif. This disrupts the ubiquitination modification process and prevents it from recognizing viral proteins, thereby preventing the initiation of innate immune response signals [32]. FMDV also negatively regulates the activation of the INF-β signaling pathway by inhibiting the phosphorylation of IRF3 and nuclear translocation through the VP1 protein [33].
Host circRNAs expression and potential role during FMDV infection are unknown. Pigs are the primary natural reservoir of FMDV and have been used as a challenge model [34]. In this study, we used RNA-Seq technology with Illumina HiSeq platform to analyze the characteristics of differentially expressed circRNAs and the biological function of parental genes in PK-15 cells after FMDV infection. RNA-reliability Seq results was verified by qRT-PCR. The results provide new clues for understanding the interaction between FMDV and host cells from the perspective of circRNAs.

Cell culture and virus infection
Porcine kidney cells (PK-15 cells) were cultured at 37 °C in a 5% CO 2 humidified atmosphere in Dulbecco's modified Eagle's medium (DMEM, Gibco, USA) supplemented with 10% heat-inactivated fetal bovine serum (FBS, Gibco, USA) and 1% antibiotic/antimycotic solution (100U/mL penicillin, 100 μg/mL streptomycin). FMDV serotype A (FMDV A/GD/MM) was provided and stored by the OIE/National Foot and Mouth Disease Reference Laboratory (Lanzhou, Gansu, China). To analyze cell response to FMDV infection, 80% confluent PK-15 cells were inoculated with FMDV or mock infected. After 8 h of infection, both the cell supernatant and precipitate were used for further analysis. All the virus-related experiments were conducted in the Biosafety Level-3 (BSL-3) laboratory of Lanzhou Veterinary Research Institute according to the standard protocols and biosafety regulations provided by the Institutional Biosafety Committee.

RNA extraction, library construction, and RNA-sequencing
Total RNA was extracted from both FMDV-infected and mock-infected cells using TRIzol reagent (Invitrogen, USA). Total RNA was quantified and the quality assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA), and 1 μg total RNA with a RIN value above 8 was used for library preparation. Pair-End index libraries were constructed according to the manufacturer's protocol (NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® ). The libraries with human circRNAs sequencing was multiplexed and loaded onto an Illumina HiSeq instrument (Illumina, USA) according to the manufacturer's instructions. RNA-seq was performed using an Illumina HiSeq to get the raw data (raw reads). Sequencing was carried out in a 2 × 150 paired-end (PE) configuration. All sequences were processed and analyzed by the Shanghai Yuanxin Biomedical Technology Company.

circRNA sequence prediction
After the libraries were sequenced, the PE reads for each sample were mapped to a reference genome (GRCh38, Ensemble 91) using TopHat2 [35,36].

Differential expression analysis of circRNAs
DE circRNAs were identified and analyzed using DESeq software based on a negative binomial distribution [37,38]. Differentially expressed (DE) genes were detected with |log 2 (fold change) |≥ 1 and P value ≤ 0.05.

Gene Ontology and Kyoto Encyclopedia of Gene and Genome pathway analysis
Gene Ontology (GO) analysis was used to annotate the biological processes involved and the functions associated with the parental genes of the circRNAs (molecular functions, cellular components, and biological processes) [39][40][41]. Kyoto Encyclopedia of Gene and Genome (KEGG) enrichment using hypergeometric test was also conducted to predict the involvement of cellular pathways targeted by circRNAs during FMDV infection [42]. The GO and KEGG pathways identified with corrected P values ≤ 0.05 were considered significantly enriched.

Protein interaction map of parental genes of differentially expressed circRNAs
Protein function related network analysis of circR-NAs genes with significant differential expression was done using the Retrieval of Interacting Genes/Proteins (STRING v10) database and web tool [43]. The DE genes were screened with |log 2 (fold change) |≥ 2 and P value ≤ 0.05.

Differentially expressed circRNAs preparation and qRT-PCR analysis
To validate the accuracy of the sequencing results and screen the partially differentially expressed circRNAs, we designed divergence primers based on the predicted sequences of the circRNAs phenology. The reliability of the sequencing results was verified by RT-qPCR based on selected circRNAs. To remove interference from linear RNA, FMDV infected group samples were treated with Ribonuclease R (RNase R), and the reaction mixture was incubated at 37 °C for 30 min and inactivated in a water bath at 72 °C for 10 min. For the control group, RNA was not treated with RNase R. After that, cDNA was synthesized using reverse transcription kit (Takara, Dalian, China) following the manufacturer's instructions. The expression of differential circRNAs was verified using the qRT-PCR assay. qRT-PCR was performed at 20 µL reaction volume, including 10 µL SYBR Green Master Mix, 1 µL PCR primers (forward and reverse, respectively), 5 µL nuclease-free water, and 3 µL cDNA. The reaction was performed at 95 °C for 2 min, followed by 45 cycles with 95 °C for 10 s and 60 °C for 10 min. The GAPDH of pigs was used as a reference house-keeping gene. All reactions were run in triplicate. The 2 −ΔΔCt method was used to measure expression level of target circRNAs.

Statistical analysis
Statistical significance analysis was performed by the Student's t-test and P values ≤ 0.05 were considered statistically significant.

Identification and classification of circRNAs
To identify differentially expressed circRNAs in FMDVinfected PK-15 cells, FMDV-infected cells and mockinfected cells were sequenced separately using the Illumina HiSeq platform. As shown in Table 1, a total of 147084824 and 140485672 circRNA raw sequence numbers were identified in FMDV-infected and mockinfected groups, respectively. After filtering and screening, 146011792 and 139140012 high quality sequences were obtained for analysis. Candidate circRNA sequences were mapped to the corresponding genomes for comparison and the FMDV-infected and mock-infected group rates were 87.18% (133109434/152684665) and 91.11% (131008908/143792709), respectively. The Q30 values for the samples were 94.83% and 95.16%. In addition, the circRNA types were identified and characterized by  Figure 1A shows that the circRNAs in the FMDV-infected and mock-infected groups were mainly derived from exons (85.17%, 87.40%), introns (8.31%, 7.77%), and intergenic regions (6.52%, 4.83%). Circos plots showed the distribution of circRNAs on chromosomes from the FMDV-infected and mockinfected groups ( Figure 1B).

Screening of differentially expressed circRNAs
With the development of sequencing technologies, high-throughput circRNA sequencing is now a standard method to evaluate circRNAs expression and it is widely used to deduce the effects of host genes in various diseases [44]. There were 1100 differentially expressed cir-cRNAs obtained from PK-15 cells after FMDV infection, of which 425 circRNAs were upregulated and 675 were downregulated. Because of the specificity of circRNAs, SRPBM (Spliced Reads per Billion Mapping) Transsplicing reads are usually used to estimate circRNA expression. We illustrated the whole expression trend and distribution of all differentially expressed circRNAs using volcano plots ( Figure 2A) and clustering heatmaps ( Figure 2B). The results showed that FMDV infection of PK-15 cells altered the intracellular circRNAs expression profile, resulting in deregulated expression of multiple circRNAs. This suggests that circRNAs may have a biological function during the cellular response to virus infection.

Functional analysis of differentially expressed circRNAs
To further reveal the potential biological functions of FMDV-induced differentially expressed circRNAs, we performed GO and KEGG pathway functional enrichment analysis for selected differentially expressed circRNA genes. GO is a widely used bioinformatics database that includes genes and gene product information for all species [45]. First, GO analysis was performed on differentially expressed circRNAs parental genes. The results showed that differentially expressed circRNAs were mainly involved in nucleic acid metabolism, cellular metabolism, and protein synthesis and modification ( Figure 3). It has been demonstrated that circRNAs can actively participate in regulating cellular metabolic processes [46,47]. Because different genes cooperate with one another to regulate biological processes, further understanding of the potential biological mechanism of differentially expressed circRNAs may be gleaned through KEGG pathway enrichment analysis. From the figures, it is clear that most of the differentially expressed circRNAs are closely associated with multiple signaling pathways including protein ubiquitination, degradation, cell cycle regulation, RNA transport, autophagy signaling, mTOR signaling, T cell receptor signaling, and nucleotide excision repair (Figure 4). Overall, we found that differentially expressed circRNA genes may regulate multiple pathways, such as nucleic acid metabolism, cellular metabolism, signal transduction, and the immune response to regulate the cellular response to viral infection. Based on gene function analysis using GO and KEGG (Tables 2 and 3), we explored the regulatory functions of these differentially expressed circRNA genes at the translation level. We screened 64 genes with significantly differentially expressed circRNAs, 8 were significantly upregulated and 56 were downregulated (p value < 0.05 and |log 2 (FC)|> 2) using the STRING for protein interaction analysis ( Figure 5). Interestingly, the results showed that these genes were primarily involved in regulating biological processes such as the protein ubiquitination, cell cycle, and nucleic acid metabolism. In vivo, ubiquitination is a key signaling and defense mechanism for host cells to detect and respond to viral infection, inducing the initiation of immune response signaling pathways against viral invasion. Based on the protein interaction network, circRNAs may play a role in the host immune response induced by FMDV infection.

Validation of the differentially expressed circRNAs by qRT-PCR
To verify the reliability of the RNA-seq data, qRT-PCR was performed to detect circRNAs expression change in PK-15 cells infected with FMDV or not. Eight candidate circR-NAs were selected to verify their reliability by qRT-PCR. The primers were designed according to the cirRNAs position (  Figure 6). These results indicated that the circRNAs in the RNA-seq datasets were reliable.

Discussion
As new rising stars of the noncoding RNA field, circR-NAs are a fascinating class of RNAs primarily found in the cytoplasm of eukaryotic cells and are considered to be the product of a reverse-splicing reaction during transcription [48]. Circular RNAs have multiple biological functions and are involved in the regulation of gene expression, disease development, and the immune response [49][50][51]. However, to date, little is known about the characteristics and functions of circRNA expression in host cells during viral infection. With the development of sequencing technology and bioinformatics, scientists have begun to unravel the roles played by circRNAs in host cells following viral infection. It has been shown that host circRNAs change after infection with certain viruses, such as porcine epidemic diarrhea virus (PEDV) [20], avian leukemia virus (ALV) [5], and transmissible gastroenteritis virus (TGEV) [52]. The function of these differentially expressed circRNAs has been analyzed by GO enrichment, KEGG pathways, and ceRNA networks to elucidate their underlying functions and mechanisms during infection.
FMD is a transmissible disease which rapidly spreads over vast areas and causes devastating effects to the livestock industry. In this study, we infected PK-15 cells with FMDV and identified 1100 differentially expressed cir-cRNAs. Through GO and STING analysis, we found that differentially expressed circRNAs were mainly involved in biological processes, such as host cell nucleic acid metabolism and cell cycle regulation induced by FMDV infection, disruption of the cell cycle, affecting the normal growth status of cells, and oncogenesis. Studies have revealed that circPLK1 promotes breast cancer cell proliferation, migration, and invasion in breast cells [53]. Overexpression of circRNACCDC66 promoted the growth and metastasis of colon cancer cells [54]. In this study, differentially expressed circRNAs, circ12: 23921820-23923094 (KPNB1), circ6: 37967643-37975235 (VPS35), circ2: 64840343-64840811 (DDX39A), and circ7: 116401757-116406130 (DICER1) are involved in biological processes, such as DNA damage repair, nucleic acid metabolism, and transport after FMDV infection. The circRNAs circ6: 96667496-96668320 (CEP192), circ6: 85556400-85558524 (RCC1), circRNA7: 39403319-39407636 (CDC5L) are involved in the regulation of the cell cycle. The above results demonstrate that the differentially expressed circRNAs in this study may be involved in the host response to FMDV infection, but the detailed functions and regulatory mechanisms of these circRNAs need to be further explored.
In addition, according to KEGG functional enrichment, most differentially expressed circRNAs were involved in biological processes such as protein ubiquitination degradation, autophagy, mTOR signaling, and cell cycle regulation, which have a close association with the immune response. Ubiquitination is a key signaling mechanism for host cells to deal with viral infections and includes post-translational modifications and proteolytic reactions. The ubiquitin-proteasome pathway is involved in regulating a variety of cellular processes, such as antigen presentation, cell cycle regulation, apoptosis, immune    responses, inflammation, and viral infections [55]. Various ubiquitination signals play an important role in the activation of the innate immune response. It has been shown that circRNAs participate in the protein ubiquitination process after viral infection. FMDV was able to disrupt the ubiquitination of proteins through multiple pathways and evade killing by the host immune system. For example, FMDV L pro could target and remove the ubiquitin-like protein, ISG15, and specifically cleave peptide bonds at the C-terminal Gly-Gly sequence to disrupt the ubiquitination modification system and inhibit the initiation of the innate immune response [32]. Orf is a worldwide zoonotic disease caused by Orf virus (ORFV). Following ORFV infection, circRNA regulates the ubiquitination process of host cells and affect the host immune response [56]. In this study, differentially expressed circ6: 79779444-79787224 (USP48) recognized and hydrolyzed the peptide bond at the C-terminal Gly residue of the  TGT ACT GAT TCC TCA TTA TGG TTA  protein and was involved in the processing of polyubiquitin and ubiquitinated proteins. It was further speculated that circRNA6: 79779444-79787224 is involved in the FMDV L pro mediated immune escape mechanism. Recently, it has been shown that circRNAs are capable of participating in the process of virus-induced autophagy in host cells. The circRNA GATAD2A was able to promote the replication of H1N1 replication by inhibiting host cell autophagy during IAV infection [21]. Overexpression of circNF1-419 regulated cellular autophagy through the PI3K-1/Akt-AMPK-mTOR and PI3K-1/Akt-mTOR signaling pathways in astrocytes, providing a new strategy for the treatment of Alzheimer's disease [57].
Other studies have shown that viruses can promote their own replication by inhibiting cellular autophagy [58]. FMDV also mediates autophagy to evade innate immunity and promote viral replication through multiple pathways. FMDV inhibits interferon production by inducing PERK-mediated autophagy [30]. FMDV VP2 can also interact with HSPB1 (Heat shock protein family B1) and activate the EIF2S1-ATF4 pathway to inhibit the AKT-MTOR signaling pathway, which inhibits autophagy in host cells to promote viral replication [31]. Heat stress stimulation can induce the synthesis of heat shock proteins which can repair the damaged proteins and degrade the unrepairable proteins as "chaperones, " maintain the stable conformation of proteins, ensure the correct folding of nascent proteins, protect cells against stress damage, enhance the tolerance of cells, and support the normal functional metabolism of cells [59,60]. In the present study, the circ2: 135356735-135368467 parental gene, HSPA4, belongs to a member of the heat shock protein family, which was significantly downregulated. We hypothesize that circ2: 135356735-135368467 may be involved in FMDV infection-mediated autophagy to evade the host immune response. Therefore, differentially expressed circRNAs in host cells may play an important role in regulating the immune response mediated by viral infection following FMDV infection. The analysis based on the above functions showed that circRNAs play important roles in regulating organism homeostasis during FMDV infection. We further validated the RNA-Seq results by circRNAs qRT-PCR. The obtained data showed a similar expression pattern compared with RNA-seq. These results indicated that the differential expression of circRNAs seen in the RNA-seq datasets was reliable. CircRNAs are usually generated by back-splicing of pre-mRNA and exhibit strong structural stability and spatiotemporal specificity [2,3,61]. Interestingly, we found that a single gene locus could be used to generate one or more circRNAs through alternative splicing. Our results show that multiple circRNA subtypes derived from the same genes were differentially expressed after FMDV infection in PK-15 cells. These circRNA parental genes play an important role in regulating the spatiotemporal expression of circRNAs.
The regulation of miRNA target gene expression by cir-cRNAs as miRNA sponges is a classical mechanism of the ceRNAs hypothesis [62,63]. Most circRNAs act as miR-NAs sponges to regulate gene expression [7,17]. ciRS-7 functions as a miR-7 sponge and contains more than 70 miRNA binding targets that inhibit miR-7 activity and promote miR-7 target gene expression [64]. The development of metastases in pancreatic cancer was associated with ciRS-7 regulating miR-7-mediated EGFR/STAT3 signaling [65]. During ORFV infection, circRNAs act as miRNA sponges to generate circRNA-miRNA-mRNA networks that indirectly regulate gene expression following ORFV infection [56]. CircRNAs have also been shown to have essential regulatory roles in colon cancer [54], breast cancer [66], and liver cancer [67]. In addition, because of the unique mode of action of circRNAs, which are stable and enriched in cells, they may be useful as molecular markers for the diagnosis of ALV and HBV infections [5,68].
In this study, we identified and analyzed the functions of differentially expressed circRNAs in host cells after FMDV infection based on GO and KEGG functional enrichment analysis, but we did not predict the miRNA target sites and miRNA target genes of the differentially expressed circRNAs. The detailed functions and regulatory mechanisms of these differentially expressed circR-NAs in viral infection or the cell cycle will be the subject of subsequent studies.
In summary, we found that FMDV infection resulted in the differential expression of circRNAs within host cells based on GO and KEGG functional enrichment analysis. The results suggest that these circRNAs are involved in the regulation of host immune response processes. In addition, this study was the first to analyze expression profiles of differential circRNAs and the biological function of parental genes derived from FMDV-infected host cells. This provides new insights for researchers to understand the mechanism underlying FMDV-host interactions from the perspective of circRNAs.