Genome-wide analysis of circular RNAs in goat skin fibroblast cells in response to Orf virus infection

Cell culture and viral infection

Goat skin fibroblast cells (GSF cells) were obtained from the Cell Bank of the Chinese Academy of Science (Kunming, China) and maintained in Dulbecco’s modified Eagle medium (DMEM; Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Invitrogen) at 37 °C and 5% CO₂. When GSF cells seeded in 10-cm diameter dish reached ∼90% confluence, the medium was removed and cells were rinsed three times with PBS. Cells were then infected with ORFV strain JLSY (tissue culture infectious dose 50% = 10^6.2/mL), a gift from Professor Guixue Hu and Research Associate Hongze Shao, at a multiplicity of infection of one. Following incubation for 1 h at 37 °C, the virus suspension was removed and cells were cultured for a further 6 h in standard medium.

RNA extraction, library construction, and sequencing

Total RNA from both ORFV infected GSF samples (OV) and uninfected GSF samples were isolated using an Ambion mir Vana miRNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA). The quality of total RNA were analyzed by Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA), and the concentration of the total RNA was quantified using a NanoDrop 2000 (Thermo Fisher Scientific, Lafayette, CO, USA). In this study, circRNA libraries were constructed as below: Ribosomal RNA was removed from 5 µg aliquots of total RNA using an Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, CA, USA). Then, the RNA fractions were fragmented and were reverse transcribed using an mRNA-Seq Sample Preparation Kit (Illumina, San Diego, CA, USA). The cDNA libraries were then paired-end sequenced using an Illumina HiSeq 4000 platform (Lc-bio, Hangzhou, China). The raw and processed data have been deposited into the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE121725. (2) Small RNA libraries: About 1 µg total RNA of each sample was used for cDNA library construction with the TruSeq Small RNA Preparation Kit (Illumina, San Diego, CA, USA) following the manufacturer’s protocol. Then, the cDNA libraries were single-end 50 bp (SE50) sequenced with an Illumina HiSeq 2500 platform (Lc-bio, Hangzhou, China). The raw and processed data have been deposited into the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE121726.

Identification and differential expression analysis of circRNAs and mRNAs, and miRNAs

Firstly, Cutadapt (Martin, 2011) was utilized to remove reads containing undetermined bases, adaptors, and low quality bases. Then sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Bowtie2 (Langmead & Salzberg, 2012) and Tophat2 (Kim et al., 2013) were used to map reads to the Capra hircus (goat) reference genome (RefSeq assembly accession: GCF_001704415.1). The matched reads of each sample were assembled using StringTie (Pertea et al., 2015). StringTie and Ballgown (Frazee et al., 2015) were utilized to evaluate the expression levels of all transcripts by Fragments per kilobase per million reads (FPKM). The dysregulated mRNAs were selected with |log₂ (fold change) | ≥ 1 and p ≤ 0.05 by R package Ballgown (Frazee et al., 2015).

Any unmapped reads were individually mapped to the goat reference genome by TopHat-Fusion (Kim & Salzberg, 2011). Then, reads mapped to the goat genome using TopHat-Fusion were analyzed by CIRCexplorer to identify candidate circRNAs (Zhang et al., 2016; Zhang et al., 2014b). The criteria were as follows: (1) GU/AG must occur at both ends of splice sites; (2) less than two mismatches; (3) more than one back-spliced junction read in at least one sample of GSF or OV group; (4) two splice sites are no more than 100 kb apart on the genome. The expression of circRNAs was calculated by the number of reads spanning back-splicing junction and FPKM was used to normalize the expression level of circRNAs. CircRNAs with |log₂ (fold change) | ≥ 1 and p ≤ 0.05 were regarded as differentially expressed by R package-edgeR (Robinson, Mccarthy & Smyth, 2014).

For miRNA analysis, ACGT101-miR (LC Sciences, Houston, Texas, USA) was used to acquire clean reads. Then, unique sequences containing 18 to 26 nucleotides were mapped to miRBase 21.0 by BLAST search to identify known and novel miRNAs. The expression level of miRNAs was normalized based on the read counts to tags per million counts (TPM). The significance standard was p ≤ 0.05.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis

GO analysis was used to determine significantly enriched GO terms (p ≤ 0.05) by hypergeometric test (Beißbarth & Speed, 2004; Gene Ontology Consortium, 2004). KEGG pathway analysis was used to explore significantly enriched pathways (p ≤ 0.05) by hypergeometric test (Kanehisa et al., 2017; Kanehisa et al., 2012). ${\displaystyle P=1-\sum _{i=0}^{m-1}\frac{\left(\genfrac{}{}{0.0pt}{}{M}{i}\right)\left(\genfrac{}{}{0.0pt}{}{N-M}{m-i}\right)}{\left(\genfrac{}{}{0.0pt}{}{N}{n}\right)}}$ N, Total number of circRNA-hosting genes; n, The number of circRNA-hosting genes with differential expression; M, The number of circRNA-hosting genes annotated to the GO term; m, The number of circRNA-hosting genes with differential expression annotated to the GO term.

Bioinformatics analysis and ceRNA network construction

The circRNA-miRNA and miRNA-mRNA interactions were predicted using TargetScan 7.0 and miRanda software. TargetScan 7.0 predicts the targets of miRNAs based on seed region homologies (Agarwal et al., 2015) while MiRanda is mainly based on the combination of free energy generated by miRNAs binding to their target genes (Betel et al., 2008). The lower the free energy, the stronger the binding. TargetScan score percentiles ≥ 50, Miranda max free energy values <−10 and Miranda score >140 were defined as the cutoff points for targets predicti. To further explore the functional role of circRNAs, a ceRNA network was constructed using Cytoscape 3.6.0 software (Shannon et al., 2003).

Validation of miRNAs

Six differentially expressed novel miRNAs (PC-3p-8215_174, PC-5p-406_14064, PC-5p-2253_1210, PC-5p-5127_361, PC-3p-10316_124 and PC-3p-4306_468) with mean TPM >30 in GSF samples or OV samples were selected for qPCR validation. First, the total RNA was reverse- transcribed using a miRNA 1st Strand cDNA Synthesis kit (Vazyme, Nanjing, China) with specific stem-loop primers (Table S1). Then, qRT-PCR assays were performed using miRNA Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) with specific forward primers (F) and the universal reverse primer (Table S1) with an ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). U6 snRNA was used as an internal control for normalization of the expression level of these miRNAs. All experiments were conducted independently three times.

Validation of circRNAs by qRT-PCR and Sanger sequencing

Differentially expressed circRNAs with mean FPKM ≥ 10 (relatively high expression ) in OV samples or GSF samples as well as possessing more than one back-spliced read in at least two replicates of OV samples or GSF samples were selected for confirmation by qRT-PCR. Thus, six down-regulated circRNAs: circRNA998, circRNA1000, circRNA1001, circRNA1684, circRNA3127, circRNA4287 and three up-regulated circRNAs: circRNA5112, circRNA7880, circRNA8565 were obtained. First, total RNA was reversely transcribed to cDNA using a Revert Aid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). Next, qPCR assays were performed with divergent primers using the ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with 2 × SYBR qPCR Mix (Aidlab, Beijing, China). 2 − ΔΔCt method was used to calculate the relative expression level of circRNAs with goat glyceraldehyde-3-phosphate dehydrogenase (GAPDH) serving as an internal control. All experiments were conducted in triplicate. The PCR products from cDNA samples were ligated into pMD-19T (Takara, Dalian, China) for Sanger sequencing to determine the back-splicing junctions.

Moreover, genomic DNA(gDNA) was extracted from GSF cells using a Blood/Cell/Tissue Genome DNA Extraction Kit (Tiangen, Beijing, China). Both cDNA and gDNA were used as templates for PCR amplification using specific convergent and divergent primers (Table 1). The PCR products were examined using 1.5% agarose gel.

Statistical analysis

Statistical significance analysis was performed by Student’s t-test, and p ≤ 0.05 was considered statistically significant.

Properties of circRNAs in ORFV-infected and uninfected GSF cells

We first performed circRNA sequencing, using the ribosomal RNA-depleted method, of three uninfected GSF samples (GSF samples) and three ORFV-infected samples (OV samples) on the Illumina HiSeq 4000 platform. We acquired an average of 86 million and 93 million raw reads for the GSF and OV groups, respectively. Clean reads, which accounted for >98.5% of the raw reads, were obtained after the removal of the low quality raw reads. The Q30 of each sample was ∼93%. Most reads (∼85%) were linearly mapped to the goat reference genome. Among the remaining unmapped reads, approximately 1% of reads from each sample were identified as back-spliced junction reads. Furthermore, we aligned the reads unmapped to goat genome from each sample to Orf virus reference genome (Orf virus NA1/11 strain under GenBank number KF234407.1). The results indicated that a total of 125,619 reads, 165,878 reads and 157,779 reads were mapped to Orf virus genome in OV-1, OV-2, and OV-3 sample, respectively. There were nearly zero reads aligned to Orf virus genome in uninfected GSF samples (Table 2).

Finally, 9,979 and 10,844 circRNAs with more than one back-spliced read existed in at least one sample of GSF or OV group were identified, respectively (Table S2). Among these circRNAs, 4,649 were shared between the groups, while 5,330 and 6,195 circRNAs were uniquely expressed in the GSF and OV samples, respectively (Fig. 1A). Interestingly, approximately 98% of the circRNAs were exonic circRNAs (ecircRNAs), while the remaining 2% were circular intronic RNAs (ciRNA) based on their location in the goat genome (Fig. 1B). These candidate circRNAs were widely distributed from chromosome 1 to chromosome 29 (Fig. 1C). In particular, chromosomes 1, 2, 3, 10, and 11 produced more than 800 circRNAs, respectively. The lengths of the circRNAs varied from 300–1,300 bp (Fig. 1D). More than 2,000 genes produced only one circRNA isoform, while 889 genes produced two different isoforms, 554 genes produced three circRNA isoforms, and 364 genes produced four isoforms, with the number of genes decreasing with the number of isoforms produced (Fig. 1E).

Additionally, an average of 40,882 transcripts were detected in GSF samples while 40,778 transcripts were identified in three OV samples, respectively. Following miRNA sequencing, 695 mature miRNAs in GSF samples and 674 miRNAs in OV samples were acquired, respectively. The mapping overview of miRNA sequencing is provided in Table S3.

Differential expression analysis of circRNAs, mRNAs and miRNAs

Based on the filtering criteria in materials and methods, a total of 151 circRNAs were identified as differentially expressed in ORFV-infected samples compared with GSF samples. Of these, 59 circRNAs were up-regulated while 92 circRNAs were down-regulated (Figs. 2A, 2B, Table S4). There were 341 differentially expressed mRNAs with 187 up-regulated mRNAs and 154 down-regulated mRNAs. The significant differences in transcript expression between OV samples and GSF samples was presented by clustered heatmap and volcano plot (Figs. 2D, 2E, Table S5). As for miRNA, 23 miRNAs were up-regulated while 33 miRNAs were down-regulated. Among the 56 differentially expressed miRNAs, 26 miRNAs were goat-derived in miRBase 21.0 and seven miRNAs (PC-3p-8215_174, PC-5p-406_14064, PC-5p-2253_1210, PC-5p-5127_361, PC-3p-10316_124, PC-3p-19472_48, PC-3p-4306_468) were first reported (Figs. 2A, 2C, Table S6).

GO and KEGG analysis for host genes of differentially expressed circRNAs

Due to the fact that the biological functions of circRNAs may be associated with their corresponding parental transcripts, GO and KEGG analyses for host genes of dysregulated circRNAs were conducted in the study. The top 20 GO terms significantly enriched (p ≤ 0.05) in molecular function, cellular component, and biological process are presented in Fig. 3A. The top six enriched GO terms in cellular component were proteinaceous extracellular matrix, intercalated disc, voltage-gated sodium channel complex, plasma membrane, T-tubule and extracellular region. The top six enriched GO terms in molecular function were growth factor binding, voltage-gated sodium channel activity, neuropilin binding, semaphorin receptor binding, chemorepellent activity and integrin binding. Moreover, regulation of inflammatory response, epithelial structure maintenance, negative regulation of insulin secretion, positive regulation of cell migration, positive regulation of ubiquitin-protein transferase activity, regulation of ion transmembrane transport were significantly enriched in the biological process subgroup (Fig. 3A, Table S7).

Next, we conducted KEGG enrichment analysis for circRNA-hosting genes. Nine pathways: Neuroactive ligand–receptor interaction, Tight junction, Rheumatoid arthritis, Transcriptional misregulation in cancers, Focal adhesion, Vascular smooth muscle contraction, Mismatch repair, Other types of O-glycan biosynthesis and Adherens junction were significantly enriched (p ≤ 0.05). Furthermore, pathways such as Platelet activation, Fc gamma R-mediated phagocytosis, Endocytosis, Cytokine-cytokine receptor interaction, and TGF-beta signaling pathway were also enriched (p > 0.05). Fc gamma R-mediated phagocytosis and cytokine-cytokine receptor interaction were involved in host immune response to pathogen infection (Fig. 3B, Table S8).

Integrated analysis of circRNAs-miRNAs-mRNAs

CircRNAs could serve as miRNA sponges indirectly regulating gene expression. Therefore, we constructed ceRNA networks based on up-regulated circRNAs, down-regulated miRNAs and up-regulated genes or down-regulated circRNAs, up-regulated miRNAs and down-regulated genes respectively to explore the biological functions of circRNAs during ORFV infection (Fig. 4). Interestingly, we found that chi-miR-103-5p_R-7, chi-miR-26b-3p, chi-miR-92a-5p, and chi-miR-122-R_1 could bind at least four circRNAs. Novel miRNA PC_3p-10316_124 was shared by both circRNA-186 and circRNA-10457 while chi-miR-2335 and circRNA8386 showed unique binding. The validated circRNA1684 could function as sponge of chi-miR-92a-5p indirectly down-regulating eleven genes and the validated circRNA-3127 could bind chi-miR-103-5P down-regulating eight genes.

Validation of miRNAs and circRNAs

The qRT-PCR results showed all six novel miRNAs could be specifically amplified. The expression levels of PC-3p-4306_468, PC-5p-5127_361 and PC-3p-10316_124 were consistent with the results of small RNA sequencing while the remaining three miRNAs did not show significant differential expression (Fig. S1).

The relative expression levels of nine circRNAs were validated by qRT-PCR using specific divergent primers. The results indicated that circRNA1001, circRNA1684 and circRN3127 were down-regulated while circRNA7880 was up-regulated, which were consistent with the RNA-seq data (Fig. 5A). As expected, the results of the agarose gel electrophoresis demonstrated that divergent primers could only amplify circRNAs from cDNA samples, while the convergent primers amplified products from both gDNA and cDNA. (Fig. 5B). Back-splicing junctions of circRNAs from cDNA samples were further validated by Sanger sequencing (Fig. 5C).