All animal experiments in our research were performed in strict accordance with the Code of Ethics of the World Medical Association (http://ec.europa.eu/environment/chemicals/lab_animals/legislation_en.htm). Experimental protocols were approved by the Animal Ethics Committee of Shanxi Agricultural University (Shanxi, China). For this study, a total of nine healthy LW and nine healthy MS male pigs were selected from the Datong Pig Breeding Farm (Shanxi, China). All pigs were kept under the same environmental conditions. Each pig was weaned and castrated at 28 days old. Three pigs each of both breeds were slaughtered at each of the following 3 developmental stages, 1, 90, and 180 days after birth. All pigs were stunned with electricity to ameliorate the suffering of the pigs before death, followed by exsanguination using transverse incision of the neck. Eight different tissues, including heart, liver, spleen, lung, kidney, pancreas, skeletal muscle, and fat, were collected, snap-frozen in liquid nitrogen, and stored at -80°C for further use. The transcriptome of longissimus dorsi at each of the 3 developmental stages was subjected to RNA sequencing analysis.
Library preparation and sequencing
The cDNA library building process was performed as follows: Total RNA was isolated from muscle tissue samples using TRIzol® Reagent (Invitrogen, USA), according to the manufacturer's instructions. RNA quantity and purity were checked using a Nanodrop 2000 (Thermo Fisher, USA), while RNA integrity was detected using agarose gel electrophoresis, with the RNA integrity number (RIN) analyzed using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). For single library preparation, the total amount of RNA should be 5 mg, at a concentration ≥ 250 ng/μL, 1.8 ≤ OD260/280 ≤ 2.2, and a RIN ≥ 7.0. Ribosomal RNAs (rRNAs) were removed from total RNA using a Ribo-Zero Magnetic kit (Epicentre, USA), and fragmented to approximately 200 bp. A TruSeqTM Stranded Total RNA Library Prep Kit (Illumina, USA) was used to prepare the cDNA library. Fragmented RNAs were used to generate first-strand cDNA using random primers. For second-strand synthesis, dTTP was substituted with dUTP in the dNTP reagent, thereby allowing the second base of the cDNA chain to contain A/U/C/G. Following end-repair and A-tailing, 150-200 bp cDNA fragments were isolated, and double-stranded cDNA was ligated to a “Y” adaptor. Single-strand cDNA was then obtained using uracil-N-glycosylase (UNG). Next, PCR amplification was performed to enrich the cDNA libraries. Finally, paired-end sequencing was used in the present study, with sequencing performed on an Illumina HiSeq4000 sequencing platform.
Quality control for raw reads and circRNA prediction
Quality control for raw reads was determined using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) software. The Q20, Q30, and GC contents of raw reads were calculated. After discarding the reads containing an adapter, reads with over 10% poly-Ns, reads with a mass value of less than 20 nt at the end of the sequence, and sequences whose length was less than 20 bp after mass pruning, remaining clean reads were aligned to the reference pig genome (Sscrofa11.1) using Hisat2 (https://ccb.jhu.edu/software/hisat2/index.shtml).
CIRI software was selected to identify circRNAs for subsequent analysis. CIRI used the BWA-MEM algorithm to perform sequence alignment and find junction reads. Junction reads, which support signals of GT-AG and alternate pairwise shear at the shear site, serve as the basis for circRNA recognition. A dynamic programming algorithm was used to detect circRNA.
Differential expression analysis of circRNA
Due to the particularity of circRNA, it is difficult to accurately obtain all circRNA reads. Therefore, the expression level of circRNA was estimated via the number of back-spliced reads. In this study, the spliced reads per billion mapping (SRPBM) method was used to estimate the circRNA expression level. The calculation formula was as follows:
SRPBM=Spliced reads/(Total mapped reads)×109.
Furthermore, edgeR (http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) was used to analyze DEcircRNAs. In this study, the screening criteria for significantly DE circRNAs were FDR < 0.05 and |log2 Fold change (FC)| ≥ 1. According to DEcircRNA SRPBM values, DEcircRNA expression pattern clustering was performed, and the distance calculation algorithm adopted; Spearman's rank correlation coefficient was used among samples, Pearson correlation coefficient was performed among circRNA, and the method of cluster was hcluster (complete algorithm).
Enrichment analysis of DEcircRNA host genes
GO analysis of DEcircRNAs was performed using GOATOOLS (https://github.com/tanghaibao/Goatools); GO terms with a corrected P-value < 0.05 were considered significantly enriched. KOBAS (http://kobas.cbi.pku.edu.cn/home.do) was used for KEGG pathway enrichment analysis; KEGG pathways with a corrected P-value < 0.05 were considered significantly enriched.
Prediction of DEcircRNA target miRNA and regulatory network construction
TargetFinder software was used to predict the DEcircRNA miRNA binding target, as well as the DEcircRNA-miRNA target regulatory relationship. Potential loci information of DEcircRNA targeting miRNA was extracted from the results according to the criteria of P ≤ 0.05 and free energy ≤ 35. The target mRNA of miRNA (DEcircRNA targeting miRNA) was further predicted. According to the binding relationship between DEcircRNA and target miRNA, and miRNA and target mRNA, the target regulatory relationship of DEcircRNA-miRNA-mRNA was obtained. Finally, Cytoscape v.3.2.1 software [64] was used to visualize each regulatory network, with default parameters adopted.
Validation of the sequencing data
To verify the reliability of the RNA sequencing data, four DEcircRNAs (circ_0094, circ_009145, circ_0017653, and circ_0014301) were randomly selected for verification. Based on NCBI reference sequences, convergent and divergent primers were designed using Primer 5.0 software to validate the existence of these circRNAs. Primer sequences are listed in Additional file 4. All primers used in this study were synthesized by Sangon (Sangon Biotech, Shanghai, China). To confirm the circRNA junction, gDNA and cDNA were used for PCR. All PCR products were sequenced by Sangon Biotech Co., Ltd. Sequence analysis was conducted using DNASTAR software (DNASTAR 7.1, http://www.dnastar.com). For RNase R treatment, 2 μg of total RNA was incubated for 20 min at 37℃ with RNase R (Lucigen Corporation, Wisconsin, USA); this was then used to synthesize cDNA for qPCR. For the control group, the same amount of RNA was incubated for 20 min at 37℃; this was then used to synthesize cDNA. qRT-PCR was performed using a SYBR PrimeScriptTM RT-PCR Kit (Takara) on an ABI-7500 (Life Technologies) under the following conditions: pre-denaturation at 95°C for 30 s, 45 cycles of 95°C for 5 s and 60°C for 34 s, one cycle of 95°C for 15 s, 60°C for 1 min, and 95°C for 30 s. All qRT-PCR reactions for each gene were performed using three biological replicates, with three replicates per experiment. 18S rRNA was used as the internal gene. The 2–ΔΔCt method was used to calculate fold changes in circRNA expression. Statistical differences among different time points were identified using ANOVA. Data were expressed as "means ± SEM". P<0.05 indicates the difference is significant. P<0.01 represents extremely significant difference.
Validation of circ_0015885
Convergent and divergent primer amplification, RNase R digestion, and qRT-PCR were used to determine the existence of circ_0015885. The circ_0015885 expression pattern at different skeletal muscle developmental stages, as well as the expression profile in various pig tissues, were detected using qRT-PCR.