Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

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

Veterinary Research volume 53, Article number: 24 (2022) Cite this article

Abstract

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 non‐coding 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 circRNAs 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 circRNAs 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 circRNAs 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 (Lpro, 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 Lpro 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.

Materials and methods

Cell culture and virus infection

Porcine kidney cells (PK-15 cells) were cultured at 37 °C in a 5% CO2 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 |log2 (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 circRNAs 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 |log2 (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.

Results

Identification and classification of circRNAs

To identify differentially expressed circRNAs in FMDV-infected PK-15 cells, FMDV-infected cells and mock-infected 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 mock-infected 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 CIRI2 sequence alignment. 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 mock-infected groups (Figure 1B).

Table 1 Summary of Illumina Hiseq sequencing data for circRNAs after infection with FMDV along with the mock group
Full size table
Figure 1
figure 1

circRNA sequencing data. (A) Categories of circRNAs in PK-15 cells infected with FMDV and mock groups. (B) Circos plots: the outline is the chromosome coverage of all circRNAs infected with FMDV along with the mock groups, whereas the inside represents the proportion of circRNAs in the reference genome.

Full size image

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 circRNAs 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) Trans-splicing 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.

Figure 2
figure 2

Volcano plot and heatmap analysis of differentially expressed circRNAs after FMDV infection. (A) The red, green, and gray dots represent upregulated, downregulated, and not significantly regulated circRNAs respectively. The x-axis indicates the fold-changes of differentially expressed circRNAs; the greater the absolute x-value, the greater the fold change. The y-axis represents the significance of the differentially expressed circRNAs; the greater the y-value, the smaller the p value. (B) Red indicates higher expression and green indicates lower expression. Each row represents a circRNA parental gene, each column represents a sample, and the tree diagram of the gene cluster is to the left.

Full size image

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.

Figure 3
figure 3

Top 30 Gene Ontology (GO) functional classifications and enrichment analysis of differentially expressed circRNAs genes. The x-axis indicates the functional description. The y-axis indicates the enrichment significance. GO functional classification from biological process (BP), cellular component (CC), and molecular function (MF).

Full size image
Figure 4
figure 4

KEGG classification and pathway enrichment analysis of differentially expressed circRNA genes. (A) KEGG classifications of differentially expressed circRNA genes. The x-axis indicates the number of differentially expressed circRNA genes and the y-axis represents KEGG terms. (B) Top 30 pathway enrichment results for differentially expressed circRNA genes. The x-axis indicates the degree of enrichment, whereas the y-axis represents the functional descriptions of the enriched pathways.

Full size image

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 |log2(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.

Table 2 Differentially expressed upregulated circRNAs (p value < 0.05)
Full size table
Table 3 Differentially expressed downregulated circRNAs (p value < 0.05)
Full size table
Figure 5
figure 5

Protein function interaction network diagram of differentially expressed circRNA genes. These genes are involved in various biological processes including nucleic acid metabolism, protein ubiquitination modifications, and cell cycle regulation.

Full size image

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 circRNAs were selected to verify their reliability by qRT-PCR. The primers were designed according to the cirRNAs position (Table 4). The FMDV infected group was compared to uninfected controls, circ12: 23,921,820–23,923,094, circ6: 37,967,643–37,975,235 were up-regulated. The circ6: 96,667,496–96,668,320, circ2:64,840,343–64,840,811, circ7: 116,401,757–116,406,130, circ2: 135,356,735–135,368,467, circ6: 85,556,400–85,558,524, circ6: 79,779,444–79,787,224 were significantly down-regulated (Figure 6). These results indicated that the circRNAs in the RNA-seq datasets were reliable.

Table 4 Primers used in differentially expressed circRNAs validation by qRT-PCR method
Full size table
Figure 6
figure 6

qRT-PCR validation of differentially expressed circRNAs genes. PK-15 cells were infected with FMDV for 8 h followed by qRT-PCR to detect differentially expressed circRNAs. Negative control (NC) indicated no infection with FMDV. The x-axis indicates the ID of the circRNAs genes, the y-axis represents the fold change after FMDV infection, and GAPDH was used as internal control.

Full size image

Discussion

As new rising stars of the noncoding RNA field, circRNAs 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 circRNAs. 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 Lpro 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 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 Lpro 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 circRNAs as miRNA sponges is a classical mechanism of the ceRNAs hypothesis [62, 63]. Most circRNAs act as miRNAs 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 circRNAs 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.

References

  1. Hsu MT, Coca-Prados M (1979) Electron microscopic evidence for the circular form of RNA in the cytoplasm of eukaryotic cells. Nature 280:339–340

    Article  CAS  PubMed  Google Scholar 

  2. Hentze MW, Preiss T (2013) Circular RNAs: splicing’s enigma variations. EMBO J 32:923–925

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Alexander L, Ulrike Z, Markus L, Christine K, Ferdinand N, Nikolaus R (2013) Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495:333–338

    Article  CAS  PubMed  Google Scholar 

  4. Zhang Y, Zhang XO, Chen T, Xiang JF, Yin QF, Xing YH, Zhu S, Yang L, Chen LL (2013) Circular intronic long noncoding RNAs. Mol Cell 51:792–806

    Article  CAS  PubMed  Google Scholar 

  5. Zhang X, Yan Y, Lei X, Li A, Zhang H, Dai Z, Li X, Chen W, Lin W, Chen F, Ma JY, Xie HM, (2017) Circular RNA alterations are involved in resistance to avian leukosis virus subgroup-J-induced tumor formation in chickens. Oncotarget 8:34961–34970

    Article  PubMed  PubMed Central  Google Scholar 

  6. Salzman J, Gawad C, Wang PL, Lacayo N, Brown PO (2012) Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS One 7:e30733

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Danan M, Schwartz S, Edelheit S, Sorek R (2012) Transcriptome-wide discovery of circular RNAs in Archaea. Nucleic Acids Res 40:3131–3142

    Article  CAS  PubMed  Google Scholar 

  8. Cocquerelle C, Mascrez B, Hétuin D, Bailleul B (1993) Mis-splicing yields circular RNA molecules. FASEB J 7:155–160

    Article  CAS  PubMed  Google Scholar 

  9. Du WW, Yang W, Liu E, Yang Z, Dhaliwal P, Yang BB (2016) Foxo3 circular RNA retards cell cycle progression via forming ternary complexes with p21 and CDK2. Nucleic Acids Res 44:2846–2858

    Article  PubMed  PubMed Central  Google Scholar 

  10. Du WW, Yang W, Chen Y, Wu ZK, Foster FS, Yang Z, Li X, Yang BB (2017) Foxo3 circular RNA promotes cardiac senescence by modulating multiple factors associated with stress and senescence responses. Eur Heart J 38:1402–1412

    CAS  PubMed  Google Scholar 

  11. Liang WC, Wong CW, Liang PP, Shi M, Cao Y, Rao ST, Tsui SK, Waye MM, Zhang Q, Fu WM (2019) Translation of the circular RNA circβ-catenin promotes liver cancer cell growth through activation of the Wnt pathway. Genome Biol 20:84

    Article  PubMed  PubMed Central  Google Scholar 

  12. Legnini I, Di Timoteo G, Rossi F, Morlando M, Briganti F, Sthandier O, Fatica A, Santini T, Andronache A, Wade M, Pietro L, Nikolaus R, Irene B (2017) Circ-ZNF609 is a circular RNA that can be translated and functions in myogenesis. Mol Cell 66:22-37.e9

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Huang H, Wei L, Qin T, Yang N, Li Z, Xu Z (2019) Circular RNA ciRS-7 triggers the migration and invasion of esophageal squamous cell carcinoma via miR-7/KLF4 and NF-κB signals. Cancer Biol Ther 20:73–80

    Article  CAS  PubMed  Google Scholar 

  14. Kleaveland B, Shi CY, Stefano J, Bartel DP (2018) A network of noncoding regulatory RNAs acts in the mammalian brain. Cell 174:350-362.e17

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Huang G, Zhu H, Shi Y, Wu W, Cai H, Chen X (2015) cir-ITCH plays an inhibitory role in colorectal cancer by regulating the Wnt/β-catenin pathway. PLoS One 10:e0131225

    Article  PubMed  PubMed Central  Google Scholar 

  16. Lukiw WJ (2013) Circular RNA (circRNA) in Alzheimer’s disease (AD). Front Genet 4:307

    Article  PubMed  PubMed Central  Google Scholar 

  17. Bachmayr-Heyda A, Reiner AT, Auer K, Sukhbaatar N, Aust S, Bachleitner-Hofmann T, Mesteri I, Grunt TW, Zeillinger R, Pils D (2015) Correlation of circular RNA abundance with proliferation–exemplified with colorectal and ovarian cancer, idiopathic lung fibrosis, and normal human tissues. Sci Rep 5:8057

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Cao S, Wei D, Li X, Zhou J, Li W, Qian Y, Wang Z, Li G, Pan X, Lei D (2017) Novel circular RNA expression profiles reflect progression of patients with hypopharyngeal squamous cell carcinoma. Oncotarget 8:45367–45379

    Article  PubMed  PubMed Central  Google Scholar 

  19. Li Y, Zheng Q, Bao C, Li S, Guo W, Zhao J, Chen D, Gu J, He X, Huang S (2015) Circular RNA is enriched and stable in exosomes: a promising biomarker for cancer diagnosis. Cell Res 25:981–984

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Chen J, Wang H, Jin L, Wang L, Huang X, Chen W, Yan M, Liu G (2019) Profile analysis of circRNAs induced by porcine endemic diarrhea virus infection in porcine intestinal epithelial cells. Virology 527:169–179

    Article  CAS  PubMed  Google Scholar 

  21. Yu T, Ding Y, Zhang Y, Liu Y, Li Y, Lei J, Zhou J, Song S, Hu B (2019) Circular RNA GATAD2A promotes H1N1 replication through inhibiting autophagy. Vet Microbiol 231:238–245

    Article  CAS  PubMed  Google Scholar 

  22. Chen TC, Tallo-Parra M, Cao QM, Kadener S, Böttcher R, Pérez-Vilaró G, Boonchuen P, Somboonwiwat K, Díez J, Sarnow P (2020) Host-derived circular RNAs display proviral activities in Hepatitis C virus-infected cells. PLoS Pathog 16:e1008346

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Liu CX, Li X, Nan F, Jiang S, Gao X, Guo SK, Xue W, Cui Y, Dong K, Ding H, Qu B, Zhou Z, Shen N, Yang L, Chen L (2019) Structure and degradation of circular RNAs regulate PKR activation in innate immunity. Cell 177:865-880.e821

    Article  CAS  PubMed  Google Scholar 

  24. Grubman MJ, Baxt B (2004) Foot-and-mouth disease. Clin Microbiol Rev 17:465–493

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Li D, Zhang J, Yang W, He Y, Ru Y, Fu S, Li L, Liu X, Zheng H (2019) Poly (rC) binding protein 2 interacts with VP0 and increases the replication of the foot-and-mouth disease virus. Cell Death Dis 10:516

    Article  PubMed  PubMed Central  Google Scholar 

  26. Jamal SM, Belsham GJ (2013) Foot-and-mouth disease: past, present and future. Vet Res 44:116

    Article  PubMed  PubMed Central  Google Scholar 

  27. Mason PW, Grubman MJ, Baxt B (2003) Molecular basis of pathogenesis of FMDV. Virus Res 91:9–32

    Article  CAS  PubMed  Google Scholar 

  28. Yang B, Zhang X, Zhang D, Hou J, Xu G, Sheng C, Choudhury SM, Zhu Z, Li D, Zhang K, Zheng H (2020) Molecular mechanisms of immune escape for foot-and-mouth disease virus. Pathogens 9:729

    Article  CAS  PubMed Central  Google Scholar 

  29. Paton DJ, Reeve R, Capozzo AV, Ludi A (2019) Estimating the protection afforded by foot-and-mouth disease vaccines in the laboratory. Vaccine 37:5515–5524

    Article  CAS  PubMed  Google Scholar 

  30. Ranjitha HB, Ammanathan V, Guleria N, Hosamani M, Sreenivasa BP, Dhanesh VV, Santhoshkumar R, Sagar BKC, Mishra BP, Singh RK, Aniket S, Ravi M, Suresh HB (2020) Foot-and-mouth disease virus induces PERK-mediated autophagy to suppress the antiviral interferon response. J Cell Sci 134:jcs240622

    Article  PubMed  Google Scholar 

  31. Sun P, Zhang S, Qin X, Chang X, Cui X, Li H, Zhang S, Gao H, Wang P, Zhang Z, Luo J, Li Z (2018) Foot-and-mouth disease virus capsid protein VP2 activates the cellular EIF2S1-ATF4 pathway and induces autophagy via HSPB1. Autophagy 14:336–346

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Swatek KN, Aumayr M, Pruneda JN, Visser LJ, Berryman S, Kueck AF, Geurink PP, Ovaa H, van Kuppeveld FJM, Tuthill TJ, Skern T, Komander D (2018) Irreversible inactivation of ISG15 by a viral leader protease enables alternative infection detection strategies. Proc Natl Acad Sci USA 115:2371–2376

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Zhang W, Yang F, Zhu Z, Yang Y, Wang Z, Cao W, Dang W, Li L, Mao R, Liu Y, Tian H, Zhang K, Liu X, Ma J, Zheng H (2019) Cellular DNAJA3, a novel VP1-interacting protein, inhibits foot-and-mouth disease virus replication by inducing lysosomal degradation of VP1 and attenuating its antagonistic role in the beta interferon signaling pathway. J Virol 93:e00588-19

    Article  PubMed  PubMed Central  Google Scholar 

  34. Cubillos C, de la Torre BG, Jakab A, Clementi G, Borrás E, Bárcena J, Andreu D, Sobrino F, Blanco E (2008) Enhanced mucosal immunoglobulin A response and solid protection against foot-and-mouth disease virus challenge induced by a novel dendrimeric peptide. J Virol 82:7223–7230

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Gao Y, Wang J, Zhao F (2015) CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol 16:4

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9:357–359

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Anders S, Huber W (2010) Differential expression analysis for sequence count data. Genome Biol 11:R106

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139–140

    Article  CAS  PubMed  Google Scholar 

  39. Tang H, Wang X, Bowers JE, Ming R, Alam M, Paterson AH (2008) Unraveling ancient hexaploidy through multiply-aligned angiosperm gene maps. Genome Res 18:1944–1954

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Young MD, Wakefield MJ, Smyth GK, Oshlack A (2010) Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol 11:R14

    Article  PubMed  PubMed Central  Google Scholar 

  41. Lu J, Peatman E, Tang H, Lewis J, Liu Z (2012) Profiling of gene duplication patterns of sequenced teleost genomes: evidence for rapid lineage-specific genome expansion mediated by recent tandem duplications. BMC Genomics 13:246

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L (2011) KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 39:W316-322

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, Mering CV (2015) STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 43:D447-452

    Article  CAS  PubMed  Google Scholar 

  44. Wang H, Yang J, Yang J, Fan Z, Yang C (2016) Circular RNAs: Novel rising stars in cardiovascular disease research. Int J Cardiol 202:726–727

    Article  PubMed  Google Scholar 

  45. Huntley RP, Sawford T, Mutowo-Meullenet P, Shypitsyna A, Bonilla C, Martin MJ, O’Donovan C (2015) The GOA database: gene Ontology annotation updates for 2015. Nucleic Acids Res 43:D1057-1063

    Article  CAS  PubMed  Google Scholar 

  46. Granados-Riveron JT, Aquino-Jarquin G (2016) The complexity of the translation ability of circRNAs. Biochim Biophys Acta 1859:1245–1251

    Article  CAS  PubMed  Google Scholar 

  47. Zou F, Ding Z, Jiang J, Lu F, Xia X, Ma X (2017) Confirmation and preliminary analysis of circRNAs potentially involved in human intervertebral disc degeneration. Mol Med Rep 16:9173–9180

    Article  CAS  PubMed  Google Scholar 

  48. Sand M, Bechara FG, Gambichler T, Sand D, Bromba M, Hahn SA, Stockfleth E, Hessam S (2016) Circular RNA expression in cutaneous squamous cell carcinoma. J Dermatol Sci 83:210–218

    Article  CAS  PubMed  Google Scholar 

  49. Li X, Liu CX, Xue W, Zhang Y, Jiang S, Yin QF, Wei J, Yao RW, Yang L, Chen LL (2017) Coordinated circRNA biogenesis and function with NF90/NF110 in viral Infection. Mol Cell 67:214-227.e7

    Article  CAS  PubMed  Google Scholar 

  50. Peng L, Chen G, Zhu Z, Shen Z, Du C, Zang R, Su Y, Xie H, Li H, Xu X, Xia Y, Tang W (2019) Correction: Circular RNA ZNF609 functions as a competitive endogenous RNA to regulate AKT3 expression by sponging miR-150-5p in Hirschsprung’s disease. Oncotarget 10:3313–3314

    Article  PubMed  PubMed Central  Google Scholar 

  51. Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, Zhong G, Yu B, Hu W, Dai L, Zhu P, Chang Z, Wu Q, Zhao Y, Jia Y, Xu P, Liu H, Shan G (2015) Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol 22:256–264

    Article  PubMed  Google Scholar 

  52. Ma X, Zhao X, Zhang Z, Guo J, Guan L, Li J, Mi M, Huang Y, Tong D (2018) Differentially expressed non-coding RNAs induced by transmissible gastroenteritis virus potentially regulate inflammation and NF-κB pathway in porcine intestinal epithelial cell line. BMC Genomics 19:747

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Lin G, Wang S, Zhang X, Wang D (2020) Circular RNA circPLK1 promotes breast cancer cell proliferation, migration and invasion by regulating miR-4500/IGF1 axis. Cancer Cell Int 20:593

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Hsiao KY, Lin YC, Gupta SK, Chang N, Yen L, Sun HS, Tsai SJ (2017) Noncoding effects of circular RNA CCDC66 promote colon cancer growth and metastasis. Cancer Res 77:2339–2350

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Schulman BA, Harper JW (2009) Ubiquitin-like protein activation by E1 enzymes. the apex for downstream signalling pathways. Nat Rev Mol Cell Biol 10:319–331

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Pang F, Zhang M, Yang X, Li G, Zhu S, Nie X, Cao R, Yang X, Zhang Z, Huang H, Li B, Wang C, Du L, Wang F (2019) Genome-wide analysis of circular RNAs in goat skin fibroblast cells in response to Orf virus infection. Peer J 7:e6267

    Article  PubMed  PubMed Central  Google Scholar 

  57. Chen D, Guo Y, Qi L, Tang X, Liu Y, Yang X, Hu G, Shuai O, Yong T, Wang D, Xie Y, Yang BB, Wu Q (2019) Circular RNA NF1-419 enhances autophagy to ameliorate senile dementia by binding Dynamin-1 and Adaptor protein 2 B1 in AD-like mice. Aging 11:12002–12031

    Article  CAS  Google Scholar 

  58. Chiramel AI, Brady NR, Bartenschlager R (2013) Divergent roles of autophagy in virus infection. Cells 2:83–104

    Article  PubMed  PubMed Central  Google Scholar 

  59. Nicoll WS, Boshoff A, Ludewig MH, Hennessy F, Jung M, Blatch GL (2006) Approaches to the isolation and characterization of molecular chaperones. Protein Express Purif 46:1–15

    Article  CAS  Google Scholar 

  60. Steel GJ, Fullerton DM, Tyson JR, Stirling CJ (2004) Coordinated activation of Hsp70 chaperones. Science 303:98–101

    Article  CAS  PubMed  Google Scholar 

  61. Barrett SP, Salzman J (2016) Circular RNAs: analysis, expression and potential functions. Development 143:1838–1847

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Taulli R, Loretelli C, Pandolfi PP (2013) From pseudo-ceRNAs to circ-ceRNAs: a tale of cross-talk and competition. Nat Struct Mol Biol 20:541–543

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Thomson DW, Dinger ME (2016) Endogenous microRNA sponges: evidence and controversy. Nat Rev Genet 17:272–283

    Article  CAS  PubMed  Google Scholar 

  64. Hansen TB, Kjems J, Damgaard CK (2013) Circular RNA and miR-7 in cancer. Cancer Res 73:5609–5612

    Article  CAS  PubMed  Google Scholar 

  65. Liu L, Liu FB, Huang M, Xie K, Xie QS, Liu CH, Shen MJ, Huang Q (2019) Circular RNA ciRS-7 promotes the proliferation and metastasis of pancreatic cancer by regulating miR-7-mediated EGFR/STAT3 signaling pathway. Hepatobiliary Pancreat Dis Int 18:580–586

    Article  PubMed  Google Scholar 

  66. Lu WY (2017) Roles of the circular RNA circ-Foxo3 in breast cancer progression. Cell Cycle 16:589–590

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Xu L, Zhang M, Zheng X, Yi P, Lan C, Xu M (2017) The circular RNA ciRS-7 (Cdr1as) acts as a risk factor of hepatic microvascular invasion in hepatocellular carcinoma. J Cancer Res Clin Oncol 143:17–27

    Article  PubMed  Google Scholar 

  68. Rao X, Lai L, Li X, Wang L, Li A, Yang Q (2021) N(6) -methyladenosine modification of circular RNA circ-ARL3 facilitates Hepatitis B virus-associated hepatocellular carcinoma via sponging miR-1305. IUBMB Life 73:408–417

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors would like to thank the anonymous editors and reviewers for their valuable comments and suggestions that proved the quality of this manuscript.

Funding

This work was supported by grants from Natural Science Foundation of Gansu Province (20JR5RA582) and National Natural Science Foundation of China (31972684).

Author information

Authors and Affiliations

  1. State Key Laboratory of Veterinary Etiological Biology, National Foot-and-Mouth Disease Reference Laboratory, Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, Lanzhou, 730046, China

    JinKe Yang, Bo Yang, Yong Wang, Ting Zhang, Yu Hao, HuiMei Cui, DengShuai Zhao, XingGuo Yuan, XueHui Chen, ChaoChao Shen, WenQian Yan, HaiXue Zheng, KeShan Zhang & Xiangtao Liu

Authors
  1. JinKe Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Bo Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Ting Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yu Hao
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. HuiMei Cui
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. DengShuai Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. XingGuo Yuan
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. XueHui Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. ChaoChao Shen
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. WenQian Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. HaiXue Zheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. KeShan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. Xiangtao Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization: KSZ, HXZ and XTL. Formal analysis: JKY, BY and YW. Funding acquisition: HXZ and KSZ. Investigation: TZ, YH, HMC, WQY and DSZ. Methodology: XGY, XHC and CCS. Writing: JKY and KSZ. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to HaiXue Zheng or KeShan Zhang.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yang, J., Yang, B., Wang, Y. et al. Profiling and functional analysis of differentially expressed circular RNAs identified in foot-and-mouth disease virus infected PK-15 cells. Vet Res 53, 24 (2022). https://doi.org/10.1186/s13567-022-01037-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13567-022-01037-w

Keywords

Veterinary Research

ISSN: 1297-9716