Landscape of tissue-specific RNA Editome provides insight into co-regulated and altered gene expression in pigs (Sus-scrofa).

RNA editing generates genetic diversity in mammals by altering amino acid sequences, miRNA targeting site sequences, influencing the stability of targeted RNAs, and causing changes in gene expression. However, the extent to which RNA editing affect gene expression via modifying miRNA binding site remains unexplored. Here, we first profiled the dynamic A-to-I RNA editome across tissues of Duroc and Luchuan pigs. The RNA editing events at the miRNA binding sites were generated. The biological function of the differentially edited gene in skeletal muscle was further characterized in pig muscle-derived satellite cells. RNA editome analysis revealed a total of 171,909 A-to-I RNA editing sites (RESs), and examination of its features showed that these A-to-I editing sites were mainly located in SINE retrotransposons PRE-1/Pre0_SS element. Analysis of differentially edited sites (DESs) revealed a total of 4,552 DESs across tissues between Duroc and Luchuan pigs, and functional category enrichment analysis of differentially edited gene (DEG) sets highlighted a significant association and enrichment of tissue-developmental pathways including TGF-beta, PI3K-Akt, AMPK, and Wnt signaling pathways. Moreover, we found that RNA editing events at the miRNA binding sites in the 3'-UTR of HSPA12B mRNA could prevent the miRNA-mediated mRNA downregulation of HSPA12B in the muscle-derived satellite (MDS) cell, consistent with the results obtained from the Luchuan skeletal muscle. This study represents the most systematic attempt to characterize the significance of RNA editing in regulating gene expression, particularly in skeletal muscle, constituting a new layer of regulation to understand the genetic mechanisms behind phenotype variance in animals.Abbreviations: A-to-I: Adenosine-to-inosine; ADAR: Adenosine deaminase acting on RNA; RES: RNA editing site; DEG: Differentially edited gene; DES: Differentially edited site; FDR: False discovery rate; GO: Gene Ontology; KEGG: Kyoto Encyclopaedia of Genes and Genomes; MDS cell: musclederived satellite cell; RPKM: Reads per kilobase of exon model in a gene per million mapped reads; UTR: Untranslated coding regions.

Introduction

One of the most prevalent forms of post-transcriptional RNA modification is the conversion of adenosine (A) nucleotides to inosine (I), mediated by the ADAR family of enzymes. Among the different types of mammalian RNA editing, adenosine-to-inosine (A-to-I) editing by double-stranded RNA-specific adenosine deaminase (ADAR) enzymes is the most common type[1] found in different organisms (Figure S1). A-to-I editing within the mRNA protein-coding region can result in codon alterations leading to the substitution of amino acids in the protein product[2], while editing in non-coding sequence usually affects transcriptional processes in specific tissue/cell [3-7]. In recent years, an increasing number of RNA editing events have been observed in miRNAs, which reveals that the editing of miRNAs influences their biogenesis, triggers their degradation, or alters the collection of the mRNAs they regulate [8,9]. Interestingly, several studies have highlighted the significance of RNA editing in pigs [7,10,11] and myogenic differentiation[12], suggesting that RNA editing events have the potential to trigger genetic variation in animals [13-16].

There are significant differences in many phenotypic traits between Chinese indigenous breeds and Western lean-type pigs, such as body size, muscle mass, and backfat thickness [17,18]. Luchuan pig is a well-known native breed in southern China known for its thin skin and excellent meat. Duroc is a famous lean-type pig breed with high muscle mass and growth rate. However, there are remarkable phenotype differences in muscle mass, growth rate, back-fat thickness and intramuscular fat content between these two breeds [19,20]. Even so, the RNA editing sites contributing to the phenotype have rarely been reported in animals. It is unclear to what extent RNA editing may be involved in tissue-specific gene expression. Advances in sequencing technologies and computational methods for identifying RESs have helped to reveal the full extent of RNA editome. Still, many questions remain unknown about the independent effect of RNA editing on miRNA binding sites i.e. the impact of A-to-I miRNA editing on miRNA biogenesis and function.

To better understand how RNA editing regulates gene expression, we first performed RNA editome profiling across ten tissues based on 57 strand-specific RNA sequencing (ss-RNA-seq) and 6 whole-genome sequencing datasets (WGS) from Luchuan and Duroc pigs. The varieties of RNA editing incidents in Luchuan and Duroc pigs’ genome were characterized, and 4552 DESs located in 1708 protein-coding genes were discovered. The genes with differential RESs were enriched in GO terms and signaling pathways associated with organ development. Furthermore, we investigated RNA editing events at miRNA binding sites in the 3′-UTR of HSPA12B, which could prevent the downregulation of the HSPA12B mRNA triggered by miRNA. Our study delineates RNA editing’s importance in regulating gene expression and provides a valuable resource for animal breeding research.

Materials and Methods

Sample collection and sequencing

The datasets of ten tissues from adult Luchuan and Duroc pigs were obtained from our previous report and are available from the China National GenBank (https://db.cngb.org/) Nucleotide Sequence Archive (CNSA) under accession number CNP000115. Briefly, three Luchuan and three Duroc adult pigs were obtained from the Institute of Animal Science, Guangxi province, China, for genome sequencing. Ten tissues were collected from each animal, respectively, representing a total of 57 RNA-Seq and 6 WGS datasets collected from organs and tissues, including heart, lung, adipose, stomach, cerebellum, liver, cerebrum, skeletal muscle, small intestine, and pancreas. Total RNA was extracted from each sample using TRIzol reagent (Invitrogen), followed by rRNA depletion and DNaseI treatment (Qiagen). For the WGS datasets, the DNeasy Blood & Tissue Kit (Qiagen) was used to isolate genomic DNA from six pigs’ ear tissues (three Luchuan and three Duroc). According to the manufacturer’s instructions, a Covaris S220 ultrasonicator (Covaris) was used to sonicate qualified DNA into 350-bp fragments. The RNA-seq libraries were constructed according to Illumina’s standard operating protocols, and the WGS libraries were generated on an Illumina HiSeq 4000 platform with 125-bp paired-end reads.

Detection of RNA editing site and validation

To rigorously detect A-to-I RESs in pig tissues, the ss-RNA-seq and WGS reads were processed before mapping to obtain high-quality reads using our state-owned Perl scripts. TopHat2[21] and BWA (v0.7.17) (http://bio-bwa.sourceforge.net/) were utilized to align clean reads to the pig reference genome (Sscrofa11.1 assembly). Clean reads that uniquely mapped (q > 10) to the reference genome were used for downstream analysis using Samtools[22] and Picard software (https://broadinstitute.github.io/picard/). Using the IndelRealigner and BaseRecalibrator toolkits from Genome Analysis Toolkit (GATK) the unique reads were subjected to local realignment. The HaplotypeCaller tool (GATK) was used to detect the variants, and snpEff (v4.3 t) was used to annotate the variants [23,24]. Next, variant quality filtering was applied, which included removing heterozygous variants; removing sites in dbSNP (version 150); removing intronic sites within 4bp of splice junctions; and removing sites within homopolymers ≥ 5bp. After that, the reads were checked with BLAT[25] and A-to-G transitions were taken as potential A-to-I RESs [16,26-29]. In total, 171,909 A-to-I RESs were identified, and thirty-five regions that harboured 80 RESs were randomly selected for validation by Sanger sequencing.

Identification of differentially editing sites (DESs)

To detect the DESs, we first compared the editing ratios between Luchuan and Duroc pigs by a 1-tailed Wilcoxon test. Welch’s t-tests were used for detecting DESs. At the same time, DEGs were identified after false discovery rate correction with an adjusted P-value (FDR-corrected) < 0.05 and an absolute fold shift ≥ 2-fold or a decrease to fold change ≤ 0.5. The DEGs that reached this threshold were preferred and subjected to downstream data analysis and functional enrichment analysis, including the association of differential RNA editing with differential gene expression, clustering, data visualization, GO annotation, and pathway analyses. RNA edited gene-sets with a nominalized P-value < 0.05 and FDR < 0.05 were defined as significantly enriched.

Data-analysis resources

For the gene expression analysis, reads were aligned to the pig reference genome (Sscrofa11.1) by TopHat2[21] using default settings. Cufflinks was used to assemble individual transcripts for comparison. The expression levels of genes were quantified using StringTie[30]. The genes differentially expressed in the two breeds were analysed using DESeq2[31] and edgeR[32]. ShinyGO – a graphical gene-set enrichment tool[33] (http://bioinformatics.sdstate.edu/go/) was used for GO annotation and pathway analyses of RNA edited genes with DESs. Cytoscape (http://www.cytoscape.org/) was used for network visualization. The annotation for major repetitive elements was downloaded from RepeatMasker (http://www.repeatmasker.org/ accessed 2020). Gene expression data sets were downloaded from Ensembl portal (https://www. ensembl.org), and miRNA targeting regions were downloaded from miRBase (release 22)[34]. Miranda software (v3.3a) (https://bioweb.pasteur.fr/packages/pack@miRanda@3.3a) was used to predict miRNA binding sites using two sets of sequences flanking the 3′-UTR regions (50 bp upstream and downstream).

Vector construction

The sequence encoding DESs of HSPA12B 3′-UTR (residues 484bp) containing miR-205 or miR-181b target sites (named HSPA12B-mutant-3′-UTR (HSPA12B-MUT-UTR) or (not-edited HSPA12B 3′-UTR (HSPA12B-WT-UTR) were amplified with appropriate primers using RNA and DNA templates. The RNA expression template was generated by annealing complementary oligonucleotides to form an RNA-coding DNA sequence and ligating the annealed DNA fragment between the Sacl and Xhol sites in pmirGLO dual-luciferase miRNA target expression vector obtained from Promega Corporation, Madison, USA. The DNA expression template was generated by annealing the DNA sequence and ligating the annealed DNA fragment between the Sacl and Xhol sites in pmirGLO dual-luciferase miRNA target expression vector obtained from Promega Corporation, Madison, USA. Both RNA and DNA expression templates were independently transformed into DH5α competent cells (Tiangen), and individual colonies were confirmed by sequencing.

Cell culture, transfection, and dual-luciferase reporter assay

To verify the target relationship between miR-205/miR-181b and DES in the 3′-UTR of HSPA12B, human embryonic kidney 293 T (HEK293T) cells were maintained in Dulbecco’s modified eagle medium (DMEM) supplemented with 10% Foetal Bovine Serum (FBS), 1% penicillin-streptomycin (PS) solution (Thermo Fisher Scientific) in a humidified incubator at 37°C and 5% CO2. HEK293T cells at 60–70% confluence were transfected with HSPA12B-WT-UTR or HSPA12B-MUT-UTR together with miR-205/miR-181b (mimics). After 48 h post-transfection, the luciferase activity was measured using Dual-Luciferase® Reporter (DLR™) Assay System (Promega, Madison, WI, USA).

HSPA12B siRNA and overexpression plasmids

In order to study the effect of HSPA12B on MDS cells, MDS cells were grown under appropriate conditions[35]. Briefly, cells were grown in DMEM supplemented with 10% FBS, 1% PS solution (Thermo Fisher Scientific) in a humidified incubator at 37°C and 5% CO2. The MDS cells at 70–90% confluence were transfected with 20 nM synthetically designed siRNA for HSPA12B (si-HSPA12B) or siRNA-negative control (si-NC) (GenePharma, Shanghai Hi-Tech, China), HSPA12B overexpression plasmid (HSPA12B-OE) or its negative control (empty vector) using LipofectamineTM 3000 reagent (Thermo Fisher Scientific, USA).

RNA extraction and quantitative real-time PCR (qRT-PCR)

The mRNA quantification of muscle-DEGs, including TM2D1, HSPA12B, CCDC86, and EIF2AK2 were assessed. Total RNA was extracted from muscle tissues using TRIzol™ LS Reagent (Thermo Fisher Scientific). A total of 500 ng of RNA was reverse-transcribed using Hiscript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd). The qRT-PCR reactions were performed using Applied Biosystems StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The expression levels of TM2D1, HSPA12B, CCDC86, and EIF2AK2 were detected using GAPDH as an endogenous control.

The relative expression level of HSPA12B and miR-181b/miR-205 (mimics) in the skeletal muscle samples were measured. Briefly, 500 ng of RNA isolated from Luchuan and Duroc skeletal muscles were reverse-transcribed using Hiscript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd). The qRT-PCR reactions were performed using Applied Biosystems StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Here, GAPDH was used as an endogenous control for HSPA12B and U6 was taken as an internal control for miR-181b/miR-205 (mimics).

The relative expression level of HSPA12B was further measured in MDS cells, with and without miR-181b/miR-205 (mimics), after 48 h of transfection. Briefly, total RNA was extracted from transfected cells using TRIzol™ LS Reagent (Thermo Fisher Scientific) according to the manufacturer’s protocol. A total of 500 ng of RNA was reverse-transcribed using Hiscript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd). The qRT-PCR reactions were performed using ChamQ Universal SYBR ® qPCR Master Mix (Vazyme Biotech Co., Ltd) on StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). GAPDH was used as an internal control.

Furthermore, the qRT-PCR of HSPA12B and miR-181b in MDS cells were quantified in anti-miR-181b (inhibitors), empty vector (control), blank well (no transfection), anti-miR-181b+si-HSPA12B, and miR-181b (mimics) templates after 48 h of transfection. In brief, total RNA was isolated using TRIzol™ LS Reagent (Thermo Fisher Scientific) according to the manufacturer’s protocol, and 500 ng of total RNA was reverse-transcribed using Hiscript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd). The qRT-PCR reactions were performed using ChamQ Universal SYBR ® qPCR Master Mix (Vazyme Biotech Co., Ltd) on Applied Biosystems StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Here, U6 was taken as an endogenous control for miR-181b and GAPDH was used as an internal control for HSPA12B. All primers or RNA-sequences used in this study are presented in Additional file 1: Table S1.

Western blot analysis

The MDS cells were transfected with the indicated oligonucleotides and plasmids for 48 h. Then, cells were lysed with Passive Lysis 5X Buffer (Promega Corporation, WI 53,711 USA). On a 10% SDS-PAGE gel, cell extracts were processed and blotted on polyvinylidene fluoride (PVDF) membrane (Millipore, Billerica, MA, USA). After blocking with 5% w/v non-fat milk, the membranes were incubated overnight at 4°C with primary antibodies including anti-PCNA (Affinity Biosciences, USA; dilution 1:5000), anti-HSPA12B (Affinity Biosciences, USA; dilution 1:1000), anti-GAPDH (Affinity Biosciences, USA; dilution 1:1000). Then, the membranes were further probed with HRP-conjugated secondary antibodies (Affinity Biosciences, USA; dilution 1:5000). The protein bands were viewed using Super ECL Detection Reagent (Yeasen Biotech Co., Ltd. Shanghai, China) on a Tanon‑5200 Chemiluminescent Imaging System (Tanon Science & Technology Co., Ltd. Shanghai, China) and the band intensity of protein was measured by ImageJ software[36].

CCK-8 assay

The MDS cells were seeded in 96-well plates at a density of 0.6 × 104 cells/ml and grown in DMEM supplemented with 10% FBS and 1% PS solution (Thermo Fisher Scientific) in a humidified incubator at 37°C with 5% CO2 for 12 h. After transfection with si-HSPA12B or si-NC, and HSPA12B-OE or empty vector (control) using (20 μM; Lipofectamine 3000 reagent (Invitrogen, Shanghai, China), cell proliferation was monitored with the CCK-8 Cell Counting Kit (Beyotime, Beijing, China)) following the manufacturer’s protocol. The absorbance at 450 nm was measured with a spectrophotometer (Tecanspark, Switzerland) after being transfected for 24, 48, 72, and 96 h.

EdU assay

Pig MDS cells were seeded in triplicate in 24-well plates (Thermo Scientific™ Nunc™) cell culture and were cultured in DMEM with Gibco™ F.B.S. for 12 h in a humidified incubator at 37°C with 5% CO2. Four groups of cells; (3 wells per group) were transfected with si-HSPA12B or si-NC, and HSPA12B-OE or empty vector (control), and another five groups of cells; (3 wells per group) were transfected with (anti-miR-181b, empty vector (control), blank well (no transfection), anti-miR-181b+si-HSPA12B, and miR-181b (mimics) using 20 μM Lipofectamine 3000 reagent (Invitrogen, Shanghai, China) for 36 h. Next, 10 mM EdU labelling medium was added to the cell culture, which was then incubated for 2 h at 37°C under 5% CO2 in a humidified incubator according to the manufacturer’s protocols (BeyoClick™ EdU-555 Cell Proliferation Detection Kit (Beyotime Biotechnology, Shanghai, China). After labelling the cells with EdU labelling medium, 4% paraformaldehyde (pH 7.4) was used at room temperature for 15 min to fix the labelled cells. The cells were washed thrice with PBS and incubated with 0.3% Triton X-100 in PBS at room temperature for 15 mins. Next, the EdU detection system was prepared according to the EdU kit manufacturer’s protocol, and the cells were incubated in the dark at room temperature for 30 mins. After a deep wash, cells were stained using 1 µg/ml of DAPI for 10 mins (Beyotime Biotechnology, Shanghai, China) followed by observation under an IX73 – inverted microscope system for live-cell imaging (Olympus, Japan). Image-Pro Plus 6.0 software (IPP 6.0) was used to estimate the differences in EdU positive cells identified by BeyoClick™ EdU-555 fluorescence.

RNA Fluorescence in situ hybridization (RNA-FISH)

Commercial FISH Kit purchased from Shanghai Gefan Biotechnology Co., Ltd (Shanghai, China) was used for the RNA-FISH assays according to manufacturer’s procedures. Briefly, pig MDS cells were fixed with 4% paraformaldehyde overnight at room temperature and washed with DEPC H2O twice for 5 mins. Cells were incubated with Proteinase K solution at 37°C for 20 mins. Then, cells were washed with PBS twice for 1 min and fixed with 4% fixation solution for 10 min at room temperature. The fixed cells were washed with PBS thrice and covered with acetic anhydride solution for 5 min at room temperature. Next, cells were washed with PBS five times and covered with pre-hybridization solution for 1 h at 65°C, followed by incubating the cells in hybridization buffer with specific probes for HSPA12B (5′-Cy3-CGCCUGCUGCCAGGAGCCCGGCUGCUCAGGUCCAGGAGCUGGUGC-3′) as well as miR-181b (5′-FAM-AACCCACCGACAGCAAUGAAUGUU-3′) for 18 h at 65°C. After a deep wash in SSC buffer, 1 µg/ml of DAPI (Beyotime Biotechnology, Shanghai, China) was used for nuclear staining. Fluorescent images were taken with an Olympus IX73 inverted microscope (Olympus, Tokyo, Japan).

Statistical analysis

Methods of 2−ΔΔCt [37,38] were used for analysing qRT-PCR results. Enrichment FDR P-value (P < 0.05) was utilized to determine enriched GO terms and KEGG pathways of genes with DESs. Also, a difference with P-value < 0.05 was taken as significant for other analyses. All data were presented as mean ± SEM or mean ± SD. Differences between the two groups were determined by the Mann-Whitney U test or a Student’s t-test. Pearson correlation coefficient was employed to determine the relationship between HSPA12B and miR-181b. P-values < 0.05 were considered statistically significant.

Results

RNA editome profiling

In this study, we performed RNA editome profiling across ten tissues using 57 ss-RNA-seq and six WGS datasets from Luchuan and Duroc to seek tissue-associated A-to-I RNA editing sites (RESs) (Fig.1Aand Additional file 2: Figure S2). A total of 177,600 RNA editing variants were identified in the pig tissue transcriptome (see Materials and Methods). As expected, A-to-G transitions were the most dominant type of RNA editing (171,909, 96.79%), followed by T-to-C variants (3,576, 2.01%, Fig. 1B). Of the 171,909 potential A-to-I RESs detected, 89,089 and 82,820 sites were observed on the forward and reverse strands, respectively (Additional file 1: Table S2 and Fig.1C). Most of these A-to-I RESs were located in major repetitive element families and types, especially in SINE/tRNA and Pre0_SS(Fig.1D and 1E).These results were consistent with previous studies in mammals [16,39]. The majority of RESs, (95.72%, 164,546/171,909) were located in introns and intergenic regions (Fig. 1F). We further compared RESs detected in this study with the Database of RNA editing in pigs (DREP); http://www.rnanet.org/editing/home.html). The DREP is a convenient online platform for free access to available RNA editing information on skeletal muscle development. Notably, 79, 159 RESs were found in the DREP (Fig 1G). This result indicated that these RESs might share common regulatory role in skeletal muscle development. To experimentally validate our analysis, we randomly picked 80 RESs from different RNA regions and amplified the genomic DNA regions for Sanger sequencing (Additional file 1: Table S3, Fig. 1H). The validation rates were 75% for editing sites in introns, 67% in missense, and 97% in 3′-UTRs (Fig. 1H). These results suggested that our RNA editome profiling was reliable for further studies.

Figure 1.

Overview of RNA editome profiling. (A) Flowchart for the identification of RESs across tissues from Duroc and Luchuan. The analyses are shown in red and green highlighted boxes. (B) Number of RNA editing variant types in the pig tissue transcriptome. (C) Total A-to-I RESs found in forward and reverse strands, respectively. (D) Distribution of RESs in major repetitive element families. (E) Distribution of RESs across the repetitive element types. (F) Distribution of RESs across different genomic regions. The numbers above the bars are the numbers of RESs. (G) The Venn diagram shows overlaps between RESs in this study and DREP database. (H) Validation of 80 RESs by Sanger sequencing and multiple alignments of reads. The validated RESs are represented by blue bars and genomic DNA mutations are represented by pink bars. The numbers inside the bars are the number of sites. In the Additional file: Table S3, validation results for these 80 RESs are reported

Differentially edited genes (DEGs) between Luchuan and Duroc pigs

Between Duroc and Luchuan pigs, 4552 differentially edited sites (DESs) were identified. These 4,552 DESs includes 649 sites in the intergenic regions and 3,903 sites in 1,708 protein-coding genes (FDR P < 0.05) (Fig. 2A and Additional file 1: Table S4). The majority of these DESs were from the small intestine (960), adipose (773), lung (747), and stomach (622), whereas only 124 DESs were detected in the skeletal muscle (Fig. 2B). Notably, the overall editing ratio in skeletal muscle was significantly lower than in other tissues (Fig. 2C), consistent with the results from a previous study in humans[40]. Also, there was a significant correlation between the editing ratio and expression levels of RNA-edited gene sets (p ≤ 0.05, Additional file 1: Table S5), suggesting that RESs potentially regulate gene expression.

Figure 2.

Differentially editing sites (DESs) between Luchuan and Duroc pigs. (A) Distribution of DESs in 3′-UTRs, 5′-UTRs, introns, missense, synonymous, and intergenic regions. Numbers above the bars are the numbers of the sites. See Additional file 1: Table S4 for all DESs identified in the present study. (B) Distribution of DESs across tissue types. The numbers above the bars are the number of sites. (C) Overall editing level across tissue types between Duroc and Luchuan. (D-E) The GO terms and Kyoto Encyclopaedia of Genes and Genomes pathways of all expressed genes harbouring DESs using ShinyGO v0.61 [33]. See Additional file 2: Table S6 and Table S7 for p-values of functional category enrichment analyses

To explore the biological functions of genes with DESs, we performed enrichment analyses of RNA-edited gene-sets for gene ontology (GO) terms and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways using ShinyGO v0.61[33]. The results revealed that genes with DESs are significantly associated with GO terms for embryo development, cell-cell signalling by Wnt, system development, cell proliferation, and cell differentiation (enrichment FDR, P < 0.05, Fig. 2D and Additional file 1: Table S6). These genes are also significantly enriched in autophagy, apoptosis, TGF-beta, PI3K-Akt, AMPK, and Wnt signalling pathways (Fig. 2E and Additional file 1: Table S7). The Wnt signaling pathway contains 19 genes, such as large tumour suppressor kinase 1 and 2 (LATS1 and LATS2), known as potential biomarkers required for tissue development [41-43]. In the TGF-beta signaling pathway, SMAD6, SMAD4, and SMAD2 genes are also associated with development [44,45]. These findings suggest that DESs in these genes may affect tissue development (Fig. 2E and Additional file 1: Table S7).

DESs in the skeletal muscle between Luchuan and Duroc pigs

In total, we found 124 DESs in the skeletal muscle between Luchuan and Duroc pigs. Of these, 65, 26, and 33 were located in introns, intergenic, and 3′-UTRs, respectively (Table 1 and Additional file 1: Table S8). As a typical regulatory strategy, miRNA can bind and regulate the stability of specific sequences in the 3′-UTR of target mRNAs [9,46,47]. Using Miranda software (v3.3a)[48], we predicted miRNA targeting regions in genes harbouring DESs in skeletal muscle (Additional file 1: Table S8). We discovered that seven of the DESs in the 3′-UTRs were located in the known miRNA targeting regions of TM2D1, HSPA12B, CCDC86, and EIF2AK2 genes (Table 2), indicating that the RESs in these genes might have a mutual effect on miRNA binding and changes in gene expression. Among the DEGs affecting miRNA binding, HSPA12B displayed up-regulated editing levels in Luchuan compared with Duroc (Fig. 3A). The gene expression changes of the RNA-edited genes, including TM2D1, HSPA12B, CCDC86, and EIF2AK2 with DESs, were examined in skeletal muscle between two breeds. The result showed that HSPA12B at mRNA level was up-regulated in Luchuan compared to Duroc (Fig. 3B and C), suggesting that RNA-editing-mutation in HSPA12B is associated with the gene expression. These findings led to the presumption that RNA editing regulates gene expression via affecting miRNA binding.

Table 1.

Genes harbouring DESs in the pig skeletal muscle

Chr.	Pos.	Strand	Repeat_type	Repeat_fam.	Gene_sym.	Anno.
1	74,651,763	+	Pre0_SS	SINE/tRNA	FOXO3	Intron
1	90,300,116	-	.	./.	SENP6	Intron
2	10,895,920	+	Pre0_SS	SINE/tRNA	CCDC86	3′-UTR
2	140,499,322	-	Pre0_SS	SINE/tRNA	ETF1	Intron
3	498,512	+	Pre0_SS	SINE/tRNA	SUN1	Intron
3	16,849,067	-	Pre0_SS	SINE/tRNA	CCT6A	Intron
3	87,276,373	+	Pre0_SS	SINE/tRNA	PSME4	Intron
3	103,090,102	+	Pre0_SS	SINE/tRNA	EIF2AK2	3′-UTR
3	103,091,518	+	.	./.	EIF2AK2	intron
3	103,093,479	+	.	./.	EIF2AK2	3′-UTR
3	103,093,488	+	.	./.	EIF2AK2	3′-UTR
3	103,093,553	+	.	./.	EIF2AK2	3′-UTR
3	103,093,561	+	.	./.	EIF2AK2	3′-UTR
4	67,478,727	+	Pre0_SS	SINE/tRNA	ARFGEF1	Intron
6	44,048,045	+	PRE1g	SINE/tRNA	GPI	Intron
6	84,856,193	-	Pre0_SS	SINE/tRNA	IFI6	Intron
6	150,653,094	+	SINE1B_SS	SINE/tRNA	TM2D1	3′-UTR
7	31,485,177	-	PRE1g	SINE/tRNA	FKBP5	Intron
8	96,972,021	+	Pre0_SS	SINE/tRNA	HSPA4L	Intron
8	120,619,275	+	Pre0_SS	SINE/tRNA	DNAJB14	3′-UTR
13	123,105,019	-	Pre0_SS	SINE/tRNA	EHHADH	3′-UTR
16	63,135,889	+	Pre0_SS	SINE/tRNA	SLU7	3′-UTR
17	31,953,686	-	Pre0_SS	SINE/tRNA	HSPA12B	3′-UTR

Chr., denotes chromosome; Pos., position; Fam., family; Sym., symbol, Anno., annotation; and UTR, untranslated region.

Table 2.

The DESs that are specifically associated with skeletal muscle in the miRNA-targeting sites

miRNA (Sus Scrofa)	Chr.	Editing site (±50bp)	Strand	Gene name	Gene ID
ssc-miR-149	Chr. 6	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-181b	Chr. 17	31,953,635–31,953,736	(-)	HSPA12B	ENSSSCG00000007147
ssc-miR-194a-3p	Chr. 2	10,895,869–10,895,970	(+)	CCDC86	ENSSSCG00000013105
ssc-miR-205	Chr. 17	31,953,635–31,953,736	(-)	HSPA12B	ENSSSCG00000007147
ssc-miR-212	Chr. 6	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-28-5p	Chr. 6	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-339	Chr. 6	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-339-5p	Chr. 6	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-345-5p	Chr. 3	103,093,502–103,093,603	(+)	EIF2AK2	ENSSSCG00000008496
ssc-miR-708-5p	Chr. 3	150,653,043–150,653,144	(+)	TM2D1	ENSSSCG00000038861
ssc-miR-7144-3p	Chr. 3	103,093,437–103,093,538	(+)	EIF2AK2	ENSSSCG00000008496

Chr., denote chromosome.

Figure 3.

The differentially edited genes (DEGs) in the pig skeletal muscle. (A) The overall editing level of DEGs located in the miRNA targeting regions. (B) The expression levels of TM2D1, HSPA12B, CCDC86, and EIF2AK2 were quantified by Reads Per Kilobase of the transcript, per Million mapped reads (RPKM). (C) The qRT-PCR results of TM2D1, HSPA12B, CCDC86, and EIF2AK2 in Luchuan and Duroc skeletal muscles. GAPDH was taken as an internal control. Significance in (A-C): *P < 0.05; **P < 0.01

RNA editing in the 3′-UTR regulate gene expression of HSPA12B via affecting miRNA binding sites for miRNA-181b and miR-205

HSPA12B is essential for cell survival and growth and can interact with multiple signaling pathways for tissue development [49,50]. The transcriptome data showed that HSPA12B was up-regulated in Luchuan when compared with Duroc pigs. Bioinformatics analysis suggested that the 3ʹ-UTR with RNA editing of HSPA12B harbours miRNAs binding sites for miR-181b and miR-205 (Table 2). The qRT-PCR results also showed that HSPA12B at mRNA level was significantly up-regulated in Luchuan pigs (Fig. 4A). Meanwhile, there was no significant difference in the expression levels of miR-181b and miR-205 between the two breeds in skeletal muscle (Fig. 4B and C), suggesting that HSPA12B abundance may not be due to miR-181b/miR-205 expression change between Luchuan and Duroc pigs. It permits hypothesizing that the DES in HSPA12B 3′-UTR might alter the targeting motif of miR-181b/miR-205 and thus, prevent miRNA-induced degradation of HSPA12B, which may result in HSPA12B up-regulation in Luchuan pigs. To validate this hypothesis, miR-181b and miR-205 were overexpressed in the skeletal muscle-derived cell line, and the expression level of HSPA12B was detected by qPCR assay. It was observed that the expression of HSPA12B at mRNA level decreased significantly as compared to the control group (Fig. 4D and E), which indicated that miR-181b and miR-205 could inhibit the expression of HSPA12B (p < 0.01 and p < 0.05, respectively).

Figure 4.

RNA editing in 3′-UTR regulates gene expression via miRNA binding sites. The relative expression level of (A), HSPA12B; (B), miR-181b; and (C), miR-205 in Luchuan and Duroc skeletal muscles. GAPDH was used as an internal control for HSPA12B and U6 for miR-181b and miR-205. Error bars indicate standard errors. (D-E) The qRT-PCR results of HSPA12B expression level in MDS cells with/without miR-181b/miR-205 transfection. GAPDH was utilized as an internal control. Error bars signify the standard errors. (F-G) The relative luciferase activity in HEK293T cells carrying edited HSPA12B-mutant-3′-UTR (HSPA12B-MUT-UTR) or not-edited HSPA12B 3′-UTR (HSPA12B-WT-UTR) in the miRNA targeting regions. Y-axis shows the relative luciferase activity. A two-tailed Wilcoxon test was used to evaluate the difference of luciferase activity using R. Error bars signify the standard errors. (H-I) Schematics of the predicted binding sites for miR-181b/miR-205 in the 3′-UTR of HSPA12B RNA editing site. The RNA mutations are shown in green letters. Significance in (A, and D-G): *P < 0.05; **P < 0.01; ***P < 0.001; ns: no significance

To further determine whether HSPA12B expression was regulated by editing the miR-181b/miR-205 binding sites, two reporter plasmids carrying a luciferase sequence with 3′-UTR of HSPA12B were generated. These plasmids had miR-181b/miR-205 targeting regions with the RNA editing site at HSPA12B. One plasmid contained a wide-type 3′-UTR editing site (HSPA12B-WT-3′-UTR), and another held mutated 3′-UTR (HSPA12B MUT −3′-UTR). These plasmids were co-transfected with miR-181b/miR-205 into the HEK293T cells, and luciferase activity was observed to determine the effect of miRNA on the expression of the reporter gene with edited or not-edited 3′-UTR of HSPA12B. It was found that the luciferase activity was downregulated to a greater degree when miR-181b and HSPA12B-WT-3′-UTR were co-transfected into the cells than with the co-transfection of HSPA12B-WT-3′-UTR (NC) (Fig. 4F). This result indicated that the luciferase activity was decreased by miR-181b for the HSPA12B with wild type 3′-UTR RNA editing site. The luciferase activity did not change when miR-181b and HSPA12B-MUT-3′-UTR were co-transfected into the cells compared to transfection of HSPA12B-MUT-3′-UTR (NC) (Fig. 4F), suggesting that the luciferase activity was not repressed for the HSPA12B with mutated 3′-UTR RNA editing site. As shown in FiG. 4G, co-transfection of miR-205 and HSPA12B-WT-3′-UTR led to a significantly decreased luciferase activity as compared with HSPA12B-WT-3′-UTR (NC). Meanwhile, the luciferase activity after co-transfection of miR-205 and HSPA12B-MUT-3′-UTR did not differ from the HSPA12B-MUT-3′-UTR (NC). These results suggest that RNA editing on the miRNA binding sites in the 3′-UTR may prevent HSPA12B mRNA degradation mediated by miR-181b/miR-205 (Fig. 4H-I).

HSPA12B promotes cell proliferation of muscle-derived satellite (MDS)

We further investigated the biological function of HSPA12B on MDS cells. Cells were cultured and stably transfected with overexpression vectors (HSPA12B-OE) and synthetically designed siRNA for HSPA12B (si-HSPA12B). Then, HSPA12B expression at mRNA and protein level was examined by qRT-PCR and western blot, respectively. We found that HSPA12B was substantially downregulated in the si-HSPA12B group compared to the si-NC and significantly up-regulated in the HSPA12B-OE group when compared to the empty vector (control) (p < 0.05; Fig. 5A–D). In si-HSPA12B and HSPA12B-OE transfected cells, there was a decrease and increase in the proliferation marker PCNA protein, respectively (Fig. 5C and D). Besides, CCK-8 and EdU assays were conducted, and the results showed that the rate of proliferation in si-HSPA12B transfected cells decreased significantly relative to HSPA12B-OE transfected cells (p < 0.05; Fig. 5E–H). These results demonstrated that overexpression of HSPA12B promotes MDS cell proliferation.

Figure 5.

The role of HSPA12B in MDS cell proliferation. (A-B) The mRNA expression level of HSPA12B in MDS cells after si-HSPA12B and HSPA12B-OE transfection. (C) The protein expression level of HSPA12B was accessed by western blot using GAPDH protein and PCNA (proliferation marker) as an internal control after si-HSPA12B transfection. Data collected from three independently replicated wells. (D) The HSPA12B protein expression level accessed by western blot using GAPDH protein as an internal control after HSPA12B-OE transfection. Data collected from three independently replicated wells. (E-F) The CCK-8 assay displayed absorbance at 450 nm in MDS cells after si-HSPA12B and HSPA12B-OE transfection. Data were collected and calculated as mean ± SD from three independently replicated wells. (G-H) The EdU positive cells in red after si-HSPA12B and HSPA12B-OE transfection. Nuclei were stained with DAPI (blue). Data were collected from three independently replicated wells. Scale bar, 100 µm. Significance in (A-B) and (E-F): **P < 0.01

HSPA12B is the target for miR-181b

Bioinformatics analyses showed that the RNA editing site in the 3′-UTR of HSPA12B contains response elements for miR-181b and miR-205. The targeting sequences of miR-181b were more abundant than miR-205 in pigs. We next chose miR-181b, which was further verified independently to be significantly up-regulated upon transfection into the MDS cells (Additional file 2: Figure S3). To detect the possible interaction between HSPA12B and miR-181b in MDS cells, we explored RNA fluorescence in situ hybridization (FISH) techniques using specific probes labelled with red and green fluorophores for the HSPA12B and miR-181b, respectively. The result showed that HSPA12B and miR-181b were co-localized in the nucleus (Fig. 6A).

Figure 6.

HSPA12B is identified as a functional target for miR-181b. (A) Co-localization of HSPA12B and miR-181b in MDS cells using RNA-FISH assay. Nuclei were stained with DAPI. HSPA12B was probed with Cy3 (red) while miR-181b was probed with FAM (green), Scale bar, 100 µm. (B) The miR-181b expression level in MDS cells after miR-181b (mimic) and anti-miR-181b (inhibitor) transfection as determined by qPCR analysis. Data are expressed as mean ± SD from three independent experiments. (C) The protein expression level of HSPA12B in MDS cells after miR-181b (mimic) and anti-miR-181b (inhibitor) transfection as determined by western blot. Data were collected from three independently replicated wells and the band intensity of HSPA12B protein expression was measured by ImageJ software. (D) HSPA12B and miR-181b relative expression levels in pig MDS cells, with each data point representing an individual sample. (E) The CCK-8 assay displayed absorbance at 450 nm in MDS cells after miR-181b (mimic) and anti-miR-181b (inhibitor) transfection. Data were collected and calculated as mean ± SD from three independently replicated wells. (F) The EdU positive cells in red after miR-181b (mimics) and anti-miR-181b transfection as determined by data collected from three independently replicated wells. Nuclei were stained with DAPI (blue). Scale bar, 100 µm. Significance in (B), (C), and (E): *p < 0.05; **P < 0.01; ***P < 0.001; ns: no significance

To further validate the relationship between miR-181b and HSPA12B, the mimics and inhibitors were used to overexpress or silence miR-181b in MDS cells, respectively (Fig. 6B). The MDS cells were cultured and transfected with mimics and inhibitors for miR-181b, and the relative protein level of HSPA12B was detected by western blot. In the MDS cell transfected with miR-181b (mimics), the relative protein level of HSPA12B decreased compared to the cell transfected with anti-miR-181b (inhibitors) (Fig. 6C). In addition, we found a significantly inverse correlation between miR-181b and HSPA12B at mRNA level in the MDS cells (y = −2.7366x+8.0561, R2 = 0.4214, p < 0.05, Pearson correlation) (Fig. 6D). This result indicates an interaction between HSPA12B and miR-181b. Furthermore, we examined the effect of miR-181b on MDS cell proliferation. Using CCK-8 and EdU assays, we observed that miR-181b knockdown by its inhibitors promotes cell proliferation while miR-181b overexpression reduced MDS cells’ proliferation (Fig. 6E and F). Altogether, these findings have shown that miR-181b regulates MDS cell proliferation through interaction with HSPA12B.

Discussion

Due to the significance of RNA editing in the regulation of gene expression [6,51], many studies on RNA editing have been conducted and documented in animal tissues [52,53] and cell lines[54]. In this study, to determine RNA editing in gene expression and/or tissue function in pigs, we performed RNA editome profiling on 57 transcriptomes across ten tissues from Duroc and Luchuan pigs. We identified 171,909 RESs, including 4,552 DESs (Fig. 1E), which can be considered as potential regulatory factors for tissue/organ development and physiological function. These 4,552 DESs include 649 sites in the intergenic regions and 3,903 sites in 1,708 protein-coding genes. Of these 3,903 sites, 587 sites were located in the 3′-UTRs, 17 in 5ʹ-UTR, 6 missense, 4 synonymous, and 3,289 in the introns (Fig. 2A and Additional file 1: Table S4). These DESs might play a role in regulating gene expression[55] reported in RNA editing studies in normal tissues[16].

DESs in introns and 3′-UTRs are enriched with known development-associated genes (Fig. 2D). In developmental contexts, however, pathways regulating tissue development were explored. For example, the hippo and Wnt signaling pathways play an essential role in tissue function [56,57]. Notably, some of the hippo-pathway targeted genes were found such as TAZ, DLG1, LAST1, LAST2, SMAD4, BMPRs, APC2, MOB1B, MOB1A, SOX5, RASSF8, RASSF1A, SMAD7, SMAD2, SMAD3, TEAD, TCF4, TCF12, TCF7L2, TCFL5, B1RC3, and B1RC6, suggesting that a comprehensive understanding of RNA modification will be of great benefit in the exploration of gene function. Also, enrichment of known developmental pathways, including TGF-beta, PI3K-Akt, AMPK, and Wnt signaling pathways, was observed (Fig. 2E). These results suggested that RESs in the 3′-UTRs and introns are involved in regulating tissue/organ development.

At the post-transcriptional level, miRNAs serve as important regulators for gene expression[58]. In animals, RNA editing at miRNA targeting sites in the 3′-UTRs is reported[8]. Also, DES in 3′-UTRs regulates skeletal muscle development via affecting miRNA binding in pigs [7,39]. In the present study, DES in the 3′-UTR of HSPA12B gene was found to regulate expression level via affecting the targeting of miR-181b in the skeletal muscle (Table 2). Previous studies have identified miR-181b in pig skeletal muscle [59,60], and miR-181b expression response to heat stress in pig[61] and other species [62,63], indicating that miR-181b is a heat-responsive molecule. Also, miR-181b is expressed in porcine embryos[64], suggesting that miR-181b may play an important role in embryonic development. In this study, RNA editing events at miR-181b/miR-205 targeting site in the 3′-UTR of HSPA12B were generated, and it was found that RNA editing on HSPA12B could prevent HSPA12B mRNA degradation mediated by miR-181b/miR-205 (Fig. 4F and G). These results suggested a possible role of RNA editing in gene expression regulation via modifying miRNA-binding sites in the 3′-UTRs.

The HSPA12B gene belongs to the heat-shock protein family and is a highly conservative protective protein found in all living organisms[65]. According to the literature, heat-shock proteins are evolutionarily conserved molecular chaperones with pivotal roles in cell survival and development [66,67]. The HSPA12B was up-regulated in the skeletal muscle of Luchuan pigs. Overexpression of HSPA12B promotes the proliferation of MDS cells (Fig. 5). These results suggest that HSPA12B-OE alone in primary cells can lead to transformation[68]. Likewise, when HSPA12B is turned off, the proliferation of MDS cells decreases. Following these observations, the overexpression of human HSPA12B in mice results in T-cell lymphomas development[69].

This study further employed RNA fluorescence in situ hybridization methods using specific probes to verify that miR-181b targets HSPA12B, probably by interacting with the 3′-UTR of HSPA12B mRNA to reduce HSPA12B expression. Besides, it was observed that the proliferation rate of MDS cells transfected with miR‐181b inhibitors markedly increased than in the cells transfected with miR-181b mimics and si-HSPA12B. This result illustrated the negative regulation of miR-181b overexpression on the proliferation of MDS cells. A potential mechanism was also observed by protein expression analysis and qRT-PCR experiment, which further demonstrated that miR-181b interacted directly with the HSPA12B protein, suggesting that the 3ʹ UTR of HSPA12B mRNA transcripts contains an important miR-181b sequence that influences the fate of HSPA12B mRNA stability and, thus, influences proteo-synthesis[70]. The mutual relationship between HSPA12B and miR-181b in MDS cells is essential for MDS cell proliferation. HSPA12B and miR-181b interactions pave the way for further study of the RNA editing phenomenon. These findings bring a more in-depth insight into the regulatory role of RNA editing in porcine MDS cells.

Conclusions

Based on the RNA/DNA sequencing technologies and bioinformatics analysis method, we have identified numerous RESs between Luchuan and Duroc pigs and found that some of these RESs overlap with those in other tissues, particularly skeletal muscles in the DREP database. Significant enrichment of RNA-edited gene-sets was found in multiple signaling pathways essential for physiological function and organ development, e.g. TGF-beta, PI3K-Akt, AMPK, and Wnt signaling pathways. More importantly, DESs in the 3ʹ-UTR of HSPA12B affect gene expression level via targeting miRNA binding sites. This exciting relationship between HSPA12B and miR-181b in the MDS cells still holds future studies. Our findings provide a basis for further investigation into the underlying role of RNA editing mutation and suggest hypotheses for further study.

Click here for additional data file.