Animals and phenotypic determination
Six barrows with 90~100 kg of live weight including three SSs and three YYs were used in this study. The SSs were purchased from the Saba pig breeding farm at Chuxiong city of the Yunnan Province, and the YYs were obtained from the pig breeding farm of the Yunnan Huijia Co. Ltd. All experimental pigs and procedures were approved by the Institutional Animal Care and Use Committee of Yunnan Agricultural University. Longissimus dorsi muscles were obtained from six pigs within fifteen minutes after slaughtered, and immediately frozen in liquid nitrogen to store for total RNA isolation. Meanwhile, longissimus dorsi tissues of 1×1×1 cm3 were processed using environment-friendly GD fixative (Servicebio, Wuhan, China) to measure muscle fiber diameter (MFD) by paraffin section. Other meat quality traits including meat color, marbling score, drip loss, water loss rate, cooked rate and intramuscular fat content (IMF) were determined according to the ‘Technical regulation for determination of pork quality (NY/T 821-2004) of Agricultural Industry Standard of China.
RNA preparation and sequencing
The total RNA was extracted from longissimus dorsi tissues using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s recommendations, and the quality of RNA was detected using 1.0% agarose gel electrophoresis and Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life, CA, USA), respectively. The RNA integrity number (RIN) was determined using the Agilent Bioanalyzer 2100 system (CA, USA). A total of 3 μg RNA per sample was used to construct small RNA sequencing library using the NEBNext Multiplex Small RNA Library Prep Set (NEB, MA, USA), and index codes were added to attribute sequences to each sample. After the clustering of the index-coded samples were produced on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia, CA, USA), and the RNA library was sequenced using Illumina Hiseq 2500 platform with 50 bp single-end (SE50) sequencing method. Similarly, Ribo-Zero RNA sequencing library was generated by 3 μg RNAs per sample. Firstly, ribosomal RNA was removed by Epicentre Ribo-Zero rRNA Removal Kit (Epicentre, WI, USA). Subsequently, the strand-specific library was generated using the rRNA-depleted RNA by NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, MA, USA) according to manufacturer’s instructions. Then, the clustering of the index-coded sample was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia, CA, USA). The library was sequenced on an Illumina Hiseq X Ten platform and 150 bp paired-end (PE150) reads were generated.
Data analysis
All raw data were filtered to generate clean data according to quality scores of all reads, and clean data were used to perform bioinformatics analysis including mapping, assembly, differentially expressed gene analysis, target gene prediction and functional enrichment analysis.
For miRNA-seq clean data, length-filtered reads were aligned to the reference genome (Sscrofa11.1, ftp://ftp.ensembl.org/pub/release-92/fasta/sus_scrofa/dna/) using Bowtie2.0 software [28]. Known miRNAs were identified using mirdeep2.0 software on the basis of miRBase database (version 22.0, http://www.mirbase.org/), and novel miRNAs were predicted according to the secondary hairpin structure using miREvo1.1 and mirdeep2.0 software [29]. The miRNA expressions were analyzed using TPM (transcripts per million reads), and differentially expressed miRNAs (DEmiRNAs) between longissimus dorsi tissues from two pig breeds were analyzed using DESeq2.0 package in Bioconductor project [30], and q value<0.05 was set as cut off criterion. DEmiRNA target genes were predicted using miRanda (http://www.microrna.org/), PITA and RNAhybrid software [31].
For Ribo-Zero RNA sequencing data, filtered paired-end reads were mapped across the pig genome (Sscrofa11.1) using Bowtie2.0 and HISAT2.0 packages with default parameters [28, 32], and the mapped reads were assembled using StringTie1.3.1 software [32]. The novel lncRNAs were identified according to the following: (1) transcripts with length shorter than 200 nt and less than 2 exons were excluded; (2) transcripts were discarded by overlapping with annotated exon region; (3) transcripts with FPKM (fragments per kilobase of transcript sequence per millions base pairs sequenced)< 0.5 were abandoned; (4) transcripts predicted with coding potential were filtered using the CNCI, CPC, Pfam-Scan and PhyloCSF tools [33-35]. Putative circRNAs were predicted using find_circ1.1 and CIRI2.1.2 packages with default parameters [36, 37], and the expressions of circRNAs were quantified using TPM. Differentially expressed lncRNAs (DElncRNAs) and circRNAs (DEcircRNAs) in longissimus dorsi tissues between SSs and YYs were screened using Ballgown and DESeq2.0 packages [30, 32], and q value< 0.05 were set as statistical significant thresholds. The targets of DElncRNA were identified based on trans (co-expression) role, and the DEcircRNA host genes were obtained according to one-to-one match between circRNA and its host genes. The functional enrichment analyses of these DEmRNAs were performed using GOseq2.12 and KOBAS2.0 software with q value< 0.05 or p value <0.05 cut off criterion [38, 39].
Competitive endogenous RNA network construction
The lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA ceRNA networks were established on the basis of ceRNA hypothesis that lncRNA and circRNA directly combine miRNA by acting as miRNA sponge to indirectly regulate the function of mRNA [21]. The lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA ceRNA networks were constructed and visualized using an open-source Cytoscape (version 3.7.0, https://cytoscape.org/) software [40].
Quantitative real-time PCR validation
To validate accuracy of RNA expression, quantitative real-time PCR (qRT-PCR) was performed for six samples including three SSs and three YYs. The cDNA was firstly reversely transcribed from total RNA using RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher , MA, USA), and qRT-PCR was implemented using ABI Step one plus platform with FastStart Universal SYBR Green Master (Roche, Basel, Switzerland) according to the manufacturer’s recommendations. Each sample was analyzed in triplicate, and the 2-ΔΔCt method was used evaluate to the RNA relative expression.