Integrated Analysis of Differentially Expressed miRNAs and mRNAs in Goat Skin Fibroblast Cells in Response to Orf Virus Infection Reveals That cfa-let-7a Regulates Thrombospondin 1 Expression

Orf is a zoonotic disease that has caused huge economic losses globally. Systematical analysis of dysregulated cellular micro RNAs (miRNAs) in response to Orf virus (ORFV) infection has not been reported. In the current study, miRNA sequencing and RNA sequencing (RNA-seq) were performed in goat skin fibroblast (GSF) cells at 0, 18, and 30 h post infection (h.p.i). We identified 140 and 221 differentially expressed (DE) miRNAs at 18 and 30 h.p.i, respectively. We also identified 729 and 3961 DE genes (DEGs) at 18 and 30 h.p.i, respectively. GO enrichment analysis indicates enrichment of apoptotic regulation, defense response to virus, immune response, and inflammatory response at both time points. DE miRNAs and DEGs with reverse expression were used to construct miRNA-gene networks. Seven DE miRNAs and seven DEGs related to “negative regulation of viral genome replication” were identified. These were validated by RT-qPCR. Cfa-let-7a, a significantly upregulated miRNA, was found to repress Thrombospondin 1 (THBS1) mRNA and protein expression by directly targeting the THBS1 3′ untranslated region. THBS1 has been reported to induce apoptosis; therefore, the cfa-let-7a-THBS1 axis may play an important role in cellular apoptosis during ORFV infection. This study provides new insights into ORFV and host cell interaction mechanisms.


Introduction
Orf, also known as contagious ecthyma, is a zoonotic disease that has led to great economic losses in livestock production globally [1]. The disease mainly affects goats and sheep, but also affects other ruminants and mammals such as musk ox, steenbok, reindeer, dog, and cat [1,2]. Moreover, people can become infected following contact with infected animals [3]. ORFV, a member of the genus Parapoxvirus, is the causative agent of Orf. ORFV has a double-stranded DNA genome of approximately 130-140 kb, encoding 132 genes [4]. The relatively conserved central regions of the ORFV genome are responsible for morphogenesis and viral replication while the highly variable terminal regions are responsible for virus virulence [5].
Various primary cell cultures and cell lines were used to isolate ORFV. Initially, primary lamb testis and primary lamb kidney cells were commonly used for ORFV isolation [6]. Primary fetal lamb muscle

Cell Culture and Viral Infection
GSF or HEK293T cells were purchased from the Cell Bank of the Chinese Academy of Science (Kunming or Shanghai, China) and maintained in Dulbecco's modified Eagle medium (DMEM; Invitrogen, Carlsbad, CA, USA), supplemented with 10% fetal bovine serum (Invitrogen). For ORFV infection, GSF cells (90% confluent) were infected with ORFV JS strain (TCID50 = 10 6.2 /mL) at a multiplicity of infection (MOI) of 1. After 1 h (h) of incubation, the supernatant was removed and cells were cultured for another 18 or 30 h.

RNA Extraction
Total RNA from uninfected GSF cells, GSF cells at 18 h.p.i and 30 h.p.i (triplicates of each group) were isolated using an Ambion mirVana™ miRNA isolation kit (Thermo Fisher Scientific, Waltham, MA, USA). The integrity and concentration of total RNA were analyzed using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and a NanoDrop™ 2000 (Thermo Fisher Scientific, Lafayette, CO, USA), respectively. Total RNA with RNA integrity number (RIN) >7 was used for high-throughput sequencing.

miRNA Sequencing and RNA-seq
As previously described, equivalent total RNA from 18 and 30 h.p.i GSF samples was used for miRNA sequencing and RNA-seq on an Illumina HiSeq 2500 or 4000 platform respectively (LC Sciences, Hangzhou, China) [29]. For miRNA analysis, raw reads from each sample were subjected to ACGT101-miR, an in-house program developed by LC Sciences, to acquire clean reads. Subsequently, unique sequences with length in 18-26 nucleotide were mapped to Capra hircus precursors in miRBase 21.0 to identify known and novel miRNAs. L/R ± n meant that the detected miRNA sequence was n base more/less than known miRNA in the left/right side. Read counts to tags per million counts (TPM) was used to normalize the expression levels of miRNAs. For RNA-seq analysis, raw data was filtered by Cutadapt to acquire clean reads [30]. Then, clean reads were aligned to the Capra hircus reference genome (Accession number: GCF_001704415.1) using the HISAT package [31]. The mapped reads of each sample were assembled using StringTie, which was then used to perform expression level for mRNAs by calculating fragments per kilobase of exon per million reads mapped (FPKM) [32].The cutoffs for DE miRNAs and DEGs were fold change ≥2 or fold change ≤0.5, and p ≤ 0.05. The raw data have been deposited in the Gene Expression Omnibus (Accession numbers: GSE141162 and GSE141163).

GO and KEGG Enrichment Analyses
To understand the biological functions and pathways of the enriched DEGs, we performed GO and KEGG pathway enrichment analyses as we previously described [29]. GO terms or pathways with p ≤ 0.05-calculated by hypergeometric test, relative to the whole genome-were significantly enriched.

Target Gene Prediction and miRNA-Gene Network Construction
MiRanda 3.3a and TargetScan 7.0 algorithms were used to predict miRNA targets in the Capra hircus genome (GCF_001704415.1). Target genes with a context score percentile of less than 50 in the TargetScan algorithm and with max free energy values > -10 in MiRanda were removed. DE miRNAs and DEGs with inverse expression were used to build miRNA-gene networks using Cytoscape 3.6.0 software [29].

RT-qPCR Validation of DEGs and DE miRNAs
For DEG validation, total RNA from each sample was used to prepare cDNA using a HiScript III first strand cDNA synthesis kit (Vazyme, Nanjing, China). Then, qPCR was performed using the ChamQ universal SYBR qPCR master mix (Vazyme) on an ABI 7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA). The miRNA first strand cDNA synthesis kit (Vazyme) and the miRNA universal SYBR qPCR master mix (Vazyme) were used for miRNA validation per manufacturer's protocols. The relative expression levels of genes or miRNAs (normalized to goat glyceraldehyde-3-phosphate dehydrogenase (GAPDH) or U6 snRNA, respectively) were calculated by the 2 −∆∆Ct method. All experiments were performed in triplicate.

Cell Transfection
GSF cells were seeded at a density of 1 × 10 5 cells/mL in a 24-well plate. Upon reaching approximately 60% confluence, cells were transfected with 100 nM cfa-let-7a_R+2 mimic or negative NC 22 control mimic (RioBio, Guangzhou, China) using lipofectamine RNAiMAX reagent (Thermo Fisher Scientific, Lafayette, Colorado, USA) per the manufacturer's protocol. After 48 h, cells were rinsed three times with PBS and lysed with an RNeasy animal RNA isolation kit (Beyotime, Shanghai, China).

Statistical Analysis
Two-tailed Students' T-test was used to evaluate the significance of the dual luciferase reporter assay and WB using GraphPad Prism 5 software. p ≤ 0.05 was considered statistically significant. Data are shown as mean ± SD from three independent experiments.

Differentially Expressed miRNAs From Intergroup Comparisons
In a previous study, we investigated changes in circular RNAs, miRNAs, and mRNAs at the early stage of ORFV infection [9]. In the current study, we focused on DE miRNAs and DEGs during the late stage of ORFV infection-18 and 30 h.p.i. Filtered raw reads yielded a total of 1465, 1438, and 1415 miRNAs in the GSF, 18, and 30 h.p.i groups, respectively ( Figure S1). Among these, 1151 miRNAs were common in all three groups. The length distribution of miRNAs in all libraries was similar, with the majority being 22 nucleotides long. At 18 and 30 h.p.i, 98 and 154 miRNAs were upregulated, respectively, while 42 and 67 miRNAs were downregulated, respectively. Compared with the 18 h.p.i group, 48 miRNAs were upregulated and 37 miRNAs were downregulated at 30 h.p.i (Figure 2A, Tables S1-S3

KEGG Enrichment Analyses of DE miRNAs From Intergroup Comparisons
To explore potential functions of miRNAs, DE miRNAs from the three comparison groups were evaluated using the TargetScan and miRanda algorithms. Predicted target genes from the analyses were then subjected to KEGG enrichment analyses. There were 121, 113, and 128 significantly enriched KEGG pathways identified in the 18 h.p.i vs. GSF, 30 h.p.i vs. GSF, and 30 h.p.i vs. 18 h.p.i comparisons, respectively (Tables S4-S6 (Table 1). After filtering low quality reads, a mean of 51,202,435, 50,443,699, and 43,420,612 clean reads were obtained from GSF, 18 h.p.i, and 30 h.p.i group, respectively. Ratios of clean reads were all above 96% in each sample. More than 93% of clean reads in GSF samples were mapped to the goat reference genome. In comparison, approximately 53% and 38% of clean reads were aligned to the goat genome in 18 and 30 h.p.i samples, respectively. The mapping percentage of clean reads to goat genome in ORFV-infected samples reduced sharply in a time-dependent manner.  Figure 4A, Tables S7-S9); Venn diagram analysis showed that 618 genes were differentially expressed in both the 18 h.p.i vs. GSF and 30 h.p.i vs. GSF comparisons and 206 genes were shared between the three intergroup comparisons ( Figure 4B).

GO Enrichment Analyses of DEGs From Intergroup Comparisons
DEGs that were enriched in the 18 h.p.i. vs. GSF comparison were mainly associated with negative regulation of apoptotic process, cell cycle, defense response to virus, immune response, and inflammatory response ( Figure 5A). DEGs that were enriched in the 30 h.p.i. vs. GSF comparison were mainly associated with cell cycle, positive regulation of apoptosis, and negative regulation of apoptosis ( Figure 5B). DEGs that were enriched in the 30 h.p.i. vs. 18 h.p.i. comparison were mainly associated with negative regulation of apoptosis, immune response, and canonical Wnt signaling pathway ( Figure 5C). The clustered heatmaps of cellular immune response genes in the 18 h.p.i vs. GSF comparison and positive and negative regulation of apoptosis genes identified in the 30 h.p.i vs. GSF comparison are shown in Figure 6.      Figure 5. In the 18 h.p.i vs. GSF comparison, cell cycle, TNF, p53, Toll-like receptor, NF-κB, and chemokine signaling pathways were significantly enriched ( Figure 7A). In the 30 h.p.i vs. GSF comparison, cell cycle, p53, TNF, phosphoinositide 3-kinase -protein kinase B (PI3K-Akt), and apoptosis pathways were significantly enriched ( Figure 7B). In the 30 h.p.i vs. 18 h.p.i comparison, PI3K-Akt, p53, cell cycle, and cytokine-cytokine receptor interaction pathways were significantly enriched ( Figure 7C).

cfa-let-7a_R+2 Target Prediction and Validation
Among the validated DE miRNAs, cfa-let-7a_R+2 was the most upregulated, therefore it was chosen for further analysis. Target gene prediction identified 55 potential targets for this miRNA. Fourteen of the identified genes had relatively high expression with FPKM >10. To verify the 14 predicted targets, cfa-let-7a_R+2 mimic or NC 22 control mimic was transfected into GSF cells. Following 48 h of incubation, RT-qPCR was conducted to detect the mRNA expression levels of the 14 targets. Results indicate that cfa-let-7a_R+2 mimic significantly reduced the mRNA expression level of beta-1,3-glucuronyltransferase 3 (B3GAT3), ring finger protein 7 (RNF7), transforming growth factor beta receptor 3 (TGFBR3), thrombospondin 1 (THBS1), and translocase of inner mitochondrial membrane 17B (TIMM17B) ( Figure 9A). Significant downregulation of RNF7, TGFBR3, THBS1, and TIMM17B was verified by RT-qPCR ( Figure 9B). Overall, THBS1 was downregulated much more than other four cfa-let-7a targets based on the above results. Furthermore, previous studies reported that THBS1 was related to apoptosis [33,34] which was an important immune mechanism used by cells against viral infection. Therefore, we first selected THBS1 for Western blot analysis. THBS1 protein expression levels were detected after 48 h following cfa-let-7a_R+2 mimic or NC 22 control mimic transfection. Results indicate that cfa-let-7a_R+2 mimic significantly reduced THBS1 protein expression compared with the control mimic ( Figure 9C,D).

Discussion
High-throughput sequencing is a powerful tool for investigating virus-host interaction. So far, this tool has been used limitedly in researching ORFV-host interactions. Therefore, in the current study, we performed miRNA sequencing and RNA-seq at 0, 18, and 30 h post ORFV infection. We identified 140 and 221 DE miRNAs at 18 and 30 h.p.i., respectively. We also identified 729 and 3961 DEGs at 18 and 30 h.p.i, respectively. GO enrichment analysis revealed that the main categories of DEG were: positive or negative regulation of apoptotic process, defense response to virus, immune response, and inflammatory response. KEGG enrichment analysis revealed association of the DEGs with TNF, p53, Toll-like receptor, NF-κB, chemokine, PI3K-Akt, and apoptosis signaling pathways.
7 DEGs related to "negative regulation of viral genome replication" were identified by RT-qPCR. Expression of ISG15 and RSAD2 displayed a sharp time-dependent increase following ORFV infection. ISG15 is an IFN-stimulated gene that encodes a ubiquitin-like protein [35,36]. ISG15 has been reported to be a broad-spectrum antiviral protein against both DNA and RNA viruses, including herpes simplex type-1, influenza A and B, HIV-1, hepatitis B and E, Ebola virus, and respiratory syncytial virus [37][38][39][40][41][42][43]. RSAD2 also known as cig5 and Viperin, is a highly conserved protein expressed in most cell types. It has been reported to be induced by double-stranded DNA, RNA, lipopolysaccharide, IFN, and a number of viruses [44]. RSAD2 has demonstrated antiviral activity against a broad range of viruses, including human cytomegalovirus, hepatitis C virus, dengue virus, influenza virus, West Nile virus, and chikungunya virus [45][46][47][48][49][50]. Whether ISG15 and RASD2 have an anti-ORFV function remains to be further studied.
Ten miRNAs DE at 30 h.p.i and 7 miRNAs DE at both 18 and 30 h.p.i were validated. One of the significantly upregulated miRNA was identified as cfa-let-7a_R+2 (hsa-let-7a_R+2). Hsa-let-7a has been reported to be downregulated in various types of human cancer including nasopharyngeal carcinoma, papillary thyroid carcinoma, lung cancer, and cervical cancer [51][52][53][54]. One study reported that hsa-let-7a inhibits migration, invasion, and tumor growth by targeting AKT2 in papillary thyroid carcinoma [51]. Another reports that hsa-let-7a inhibits proliferation and induces apoptosis by targeting EZH2 in nasopharyngeal carcinoma cells [52]. Moreover, hsa-let-7a has been shown to elevate p21WAF1 levels by targeting NIRF and suppressing the growth of A549 lung cancer cells [53].
In the current study, cfa-let 7a was identified to suppress THBS-1 mRNA and protein expression by directly targeting THBS1 3 UTR. THBS-1 is a multifunctional extra-cellular matrix glycoprotein secreted by multiple types of cells including endothelial cells, monocytes, and fibroblasts [33,55]. THBS-1 is an endogenous inhibitor of angiogenesis [56]. A previous study reported that THBS-1 inhibits angiogenesis and induces endothelial cell apoptosis [33]. Zhu et al. demonstrated that miR-222 inhibits apoptosis in porcine follicular granulosa cells by suppressing the expression of THBS1 [34]. Based on previous studies and the current results, we postulate that cfa-let-7a suppresses cellular apoptosis by targeting THBS1, which would be beneficial for ORFV replication in GSF cells. This assumption needs to be further confirmed in the future study.
Author Contributions: F.P. and F.W. conceived and designed the experiments. F.P., X.W. and Z.C. performed the experiments. F.P., X.W., Z.Z., M.Z. and C.W. prepared figures and tables. F.P., X.W., Z.C. and L.D. analyzed the data. Z.Z., M.Z., X.Y. and Q.A. contributed reagents/materials. F.P. and F.W. authored the manuscript. All authors have read and agree to the published version of the manuscript.