Genome-wide profiling of Sus scrofa circular RNAs across nine organs and three developmental stages.

The spatio-temporal expression patterns of Circular RNA (circRNA) across organs and developmental stages are critical for its function and evolution analysis. However, they remain largely unclear in mammals. Here, we comprehensively analysed circRNAs in nine organs and three skeletal muscles of Guizhou miniature pig (S. scrofa), a widely used biomedical model animal. We identified 5,934 circRNAs and analysed their molecular properties, sequence conservation, spatio-temporal expression pattern, potential function, and interaction with miRNAs. S. scrofa circRNAs show modest sequence conservation with human and mouse circRNAs, are flanked by long introns, exhibit low abundance, and are expressed dynamically in a spatio-temporally specific manner. S. scrofa circRNAs show the greatest abundance and complexity in the testis. Notably, 31% of circRNAs harbour well-conserved canonical miRNA seed matches, suggesting that some circRNAs act as miRNAs sponges. We identified 149 circRNAs potentially associated with muscle growth and found that their host genes were significantly involved in muscle development, contraction, chromatin modification, cation homeostasis, and ATP hydrolysis-coupled proton transport; moreover, this set of genes was markedly enriched in genes involved in tight junctions and the calcium signalling pathway. Finally, we constructed the first public S. scrofa circRNA database, allowing researchers to query comprehensive annotation, expression, and regulatory networks of circRNAs.

1. Introduction

Circular RNAs (circRNAs) were generally believed to be linear until their circular form was observed via electron microscopy in 1979. CircRNA is formed by exon-scrambling, in which a splice donor is spliced to an upstream acceptor rather than a downstream acceptor during transcription. While exons are most common, some introns can also form circRNAs during transcription. Usually, circRNAs are expressed at a relatively low level. CircRNAs have been identified in large numbers by deep RNA sequencing and are believed to be widely expressed in eukaryotic cells. Most circRNA is flanked by GT/AG splice sites reflecting back-splicing. CircRNA back-splicing locations are in nearly the same locations in orthologous genes and exons in animals, showing that circRNAs are conserved in animals, and plants. Although many circRNAs have been discovered in mouse and human, little is known about circRNAs in other mammals.

CircRNA expression is spatio-temporally regulated in Drosophila brains, pig brains, and mouse pre-implantation embryos. CircRNAs also show tissue-specific expression in the heart and lung during human fetal development. It is believed that circRNAs bind miRNAs and thus regulate gene expression at the post-transcriptional level. For example, CDR1as functions as an efficient miRNA sponge, binding miR-7 and allowing mRNAs to escape degradation following miRNA binding.,Sry can form a circRNA that is also proposed to act as a miRNA sponge. In addition to regulating mRNA translocation by acting as a miRNA sponge, circRNA can also regulate transcription through specific RNA–RNA interactions. A recently described subclass of circRNAs, EIciRNAs, contains introns among exons after RNA circling. In human cells, some EIciRNAs associate with RNA polymerase II, thus interacting with U1 snRNP and enhancing parental gene expression in cis. The functions of circRNAs remain largely unclear. Investigating the spatio-temporal distribution of circRNAs is critical to understanding their biological functions. However, little is known about the expression patterns of circRNAs across organs and developmental stages in mammals.

The domestic pig (Sus scrofa) has had a close and complex relationship with humans for at least 10,000 years. The pig is an important domesticated farm animal that provides protein for humans, as well as a non-rodent animal model used in biomedical research. In comparison with traditional rodent models, the pig is more similar to humans in body size, growth, development, immunity, physiology, and metabolism, as well as genomic sequence. Recent studies of the S. scrofa genome and transcriptome have increased the utility of this species as a biomedical model., However, in comparison with mice and humans, knowledge regarding the pig transcriptome across organs and developmental stages is very limited. Thousands of circRNAs have been discovered in fetal pig brains; however, the identities and expression profiles of S. scrofa circRNAs remain largely uncharacterized.

We performed genome-wide analysis of circRNA across nine organs (heart, liver, spleen, lung, kidney, ovarium, testis, muscle, and fat) and three skeletal muscles at multiple developmental stages (0, 30, and 240 days after birth) from Guizhou miniature pigs (a native Chinese pig in the Guizhou province), which are widely used as an animal model in biomedical research. We systematically identified S. scrofa circRNAs and analysed their molecular properties, sequence conservation, spatio-temporal expression patterns, potential functions, and binding sites for miRNAs. Finally, we constructed the first public S. scrofa circRNA database.

2. Materials and methods

Nine organs, including the heart, liver, spleen, lung, kidney, ovarium, testis, skeletal muscle, and fat, were collected from three male and female Guizhou miniature pigs at postnatal day 240. We also collected longissimus muscle samples from three Guizhou miniature pigs at postnatal days 0 and 30. Total RNA was isolated from each sample by TRIZOL reagent according to the manufacturer's protocols, after which total RNA was pooled by tissue type and age for RNA-seq analysis. Total RNA was reverse-transcribed according to the manufacturer’s instructions. Paired-end reads 100 bp in length were sequenced by an Illumina HiSeq 2500. RNA-seq data were deposited in the Gene Expression Omnibus (accession codes GSE73763).

Reads were first mapped to the pig genome by Bowtie2, after which the mapped reads were filtered out by SAMtools. For the unmapped reads, 20 mers were cut off from both ends of each paired-end read, after which the remaining 60 mers were aligned to the pig genome in reverse orientation by Bowtie2. We used the find_circ package, to identify circRNA. The workflow is shown in Supplementary Fig. S5.

In this study, a circRNA was designated as conserved when it was produced from orthologous genes in different species. We extracted pig-to-human and pig-to-mouse orthologous gene tables from BioMart (www.biomart.org). Human and mouse circRNAs were identified by a previous study.,

CircRNA expression was quantified by calculating the number of junction spanning reads per million reads of each sample. Junction spanning reads that were less than 5 reads supported were excluded when determining circRNA expression levels. Differentially expressed circRNAs were identified by the chi-square test with P ≤ 0.05 and had fold-change ≥ 2 (or ≤ 0.5) between any two tissues samples.

Four-hundred and eleven mature pig miRNA sequences were obtained from miRBase. We used microtar, Miranda, and RNAhybrid to identify circRNA binding sites for miRNAs. To decrease the proportion of false-positive results, intersections predicted by all three programs were treated as the final set of target miRNAs.

The circRNA database is composed of a web interface and a SQLite database engine, which is used to store and manage all data. The data processing programs are written in Python (version 3.4.4). The web interface is built by Django (version 1.8.13) a Python Web framework. Google© charts (https://developers.google.com/chart/) and D3.js (https://d3js.org/) were used for interface development. The circRNA database is coded as a Django project and is deployed by uWSGI (https://uwsgi-docs.readthedocs.io). Web services were built using nginx (version 2.10.0), an HTTP and reverse proxy server.

3. Results

3.1. Identification and characteristics of S. scrofa circRNAs

To perform comprehensive profiling of circular RNA in pigs, we carried out total RNA sequencing in nine adult organs (including fat, heart, kidney, liver, lung, spleen, testis, ovarium, and skeletal muscle) and three skeletal muscles at 0, 30, and 240 days after birth. In this study, the identified circRNAs were flanked by GT/AG splice sites and had unambiguous break points. The genomic distance between two splice sites was less than 100 kb, and each circRNA was supported by at least 2 reads in each sample. We identified 5,934 unique circRNAs in all assessed biological tissues (Supplementary Table S1).

The analysis of chromosome distribution showed that circRNAs are widely and unevenly transcribed from S. scrofa chromosomes (SSCs). More circRNAs were generated from SSC1 and SSC6 in comparison with the other SSCs. Only three circRNAs were detected on SSC-Y (Fig. 1A). The number of circRNAs was not significantly related to chromosome length or the number of genes contained in each chromosome (Supplementary Fig. S1). Most circRNAs contain multiple exons, and some circRNAs also retain introns, as indicated in a previous study. According to the pig genes and splice junction annotations in the Ensembl database (release 78, ftp://ftp.ensembl.org/pub/release-78/gtf/sus_scrofa), we found that most circularizing events typically encompass less than 5 exons (average exons = 4.8). More than 90% of circRNAs are formed by multiple exons, while only 9.9% are formed by a single exon (Fig. 1B). This finding is in agreement with observations in mice.

Figure 1

Identification and feature of circRNAs in pigs. (A) Chromosome distribution of circRNAs, the curve line represents the density of circRNAs along chromosome. (B) Distribution of exon number for circRNAs and mRNA. (C) Genomic location of circRNAs. (D) GC contents of circRNAs comparing with mRNA, miRNA and lincRNA.

CircRNAs are mostly comprised of coding sequence (CDS, 68.40%) and intergenic regions (21.93%), whereas a smaller fraction (10%) also contain 5′-untranslated regions (UTRs), 3′-UTRs, or both UTRs (Fig. 1C), in agreement with a previous report. The distribution of circRNA exonic sequence lengths resembles that of human circRNA, with a median length of 509 nucleotides, whereas that of human circRNA is 547 nucleotides. Most S. scrofa circRNAs are shorter than 2 kb (4,243/5,934 = 71.5%), similar to mouse circRNAs. The GC content of S. scrofa circRNA is similar to that of other small RNAs, with an average GC content of 42.67% (Fig. 1D).

3.2. Sequence conservation among pig, human, and mouse circRNAs

Several studies suggest that some circRNAs are evolutionarily conserved among humans and mice., Indeed, 81 of 1,903 mouse circRNAs were mapped to human circRNAs in a report by Memczak. Using the liftOver tool, another study found that 28.5% of mouse circRNAs are conserved in humans, while 33.3% of these circRNAs have no human homolog. To investigate circRNA sequence conservation between pigs and other mammals, we compared pig circRNAs to human and mouse circRNAs with clear one-to-one orthologs (http://asia.ensembl.org/biomart/martview/). We found that 20.20% of pig circRNAs have human orthologs, whereas 16.96% of pig circRNAs have mouse orthologs (Fig. 2A). These observations are comparable with previous reports.

Figure 2

Molecular characteristics of S. scrofa circRNA. (A) Sequence conservation of S. scrofa circRNA compared with human and mouse ones, respectively. The percentage means that if there is a circRNA in pig, there is also have an otology circRNA in human or mouse. (B) Sequences identity of flanking introns of a circRNA 5' end upstream and a 3' end downstream. (C) Distribution of flanking intron sequences lengths for circRNA. (D) Distribution of sequences identity of each intron pair that flanked circRNA.

To assess whether pig and human circRNAs arise from orthologous exons, we used whole-genome alignments to identify regions of the human genome that correspond to pig circRNAs and quantified the degree to which pig circRNAs overlapped with these regions. We found that 29.4% of pig circRNAs have orthologous exons in humans. In mice, 28.6% of circRNAs overlap with human circRNAs. Blast analysis suggested that 1,510 (1510/5934 = 25.45%) and 5,189 (5189/5934 = 87.44%) S. scrofa circRNAs have orthologs in mice and humans, respectively (Supplementary Tables S2 and S3). These findings indicate that pig, human, and mouse circRNAs are modestly conserved and often generated from orthologous genes. Moreover, our sequence conservation analysis suggests that circRNA might have conserved functions in pigs, mice, and humans.

3.3. S. scrofa circRNAs are flanked by long introns

The sequence conservation of circRNAs among species indicated that their mechanisms of biogenesis might also be conserved. Previous studies suggested that most circRNAs are derived from exons with long flanking introns, which usually contain reverse complementary matches (RCMs) that facilitate back-splicing of the enclosed circRNAs., To determine whether introns flank pig circRNAs, we analysed 5kb sequences length of both 5′ upstream and 3′ downstream sequences flanking pig circRNAs. We observed that the majority of circularized exons are flanked by upstream and downstream introns. We found that 60.4% of 5′ sequences upstream of circRNAs contain introns, whereas 62.4% of 3′ sequences downstream of circRNAs contain introns. Moreover, 29.1 and 28.0% of the assessed regions upstream and downstream of circRNAs are intergenic regions. Few sequences flanking circRNAs are CDS, 5′-UTRs, or 3′-UTRs (Fig. 2B). We found that 2,641 circRNAs are flanked by intron sequences on the 5′ and 3′ ends. The median length of the flanking introns on both sides of the identified circRNAs is 5,040 bp (Fig. 2C). Our analysis suggests that long flanking introns are an intrinsic determinant of RNA circularization in pigs. Mammalian circRNAs have been reported to have relatively long flanking introns., We aligned the introns of each intron pair that flanked a circRNA using the Basic Local Alignment Search Tool (BLAST). We observed that 1,810 circRNAs share RCMs of at least 100 bp in pigs (Fig. 2D). The lengths of long introns are not randomly distributed among circularized exons in the genome. Short intronic repeat sequences and complementary sequences facilitate RNA circularization., In our study, few circRNAs flanking intronic repeat sequences were detected, but the RCMs detected in our experiments are longer than previously reported RCMs; these characteristics may be typical of S. scrofa circRNAs.

3.4. circRNA isoforms and diversity of circularization pattern

Most genes generated one or two circRNAs, but some genes yielded multiple distinct circularized products (Fig. 3A). We identified 4,928 circRNAs generated from 2,538 annotated host genes. Among the identified circRNAs, 85.7% are produced by ‘hot-spot’ host genes that generate more than two circRNAs. This finding has been reported in other studies of circRNAs.,, A striking example is the chromodomain helicase DNA binding protein 2 gene (CHD2), which is involved in epigenetic regulation of chromatin structure and may generate as many as 40 distinct circRNAs (at least two unique back-spliced reads). CircRNA isoforms provide insight into the mechanism of circRNA biogenesis and the manner in which circRNAs are regulated. It is notable that, because de novo circRNA detection was performed without prior knowledge of annotated exons, 1,006 of the identified circRNAs fall outside the genomic regions of annotated genes.

Figure 3

Generation model and examples of circRNAs isoform in pig. (A) Distribution of circular RNAs among genes; (B) The genome location of three deduced donors for one acceptor to generate circRNA in DIS3L2. (C) The genome location of two deduced acceptors for one donor to generate circRNA in RBM25. (D) Expression profile of three circRNA isoforms for DIS3L2 in nine tissues. (E) ssc-ciR-03122 is abundantly expressed in kidney. (F) Tissue and age-dependent differential expression for two circRNA isoforms in PTGES3.

There are structural complexities among well-expressed circRNAs. circRNAs produced from back-splicing have splicing signals at their junctions. Introns spliced by the major spliceosome usually contain the GT dinucleotide at their 5′ end (the splice donor) and the AG dinucleotide at their 3′ end (the splice acceptor). Most circRNAs contain one or a few exons; however, many circRNA are supported by abundant back-spliced reads that traverse multiple exons. More than one GT donor site often exists for a given AG acceptor splicing site. We found that this situation is common when GT donors are more than one exon away from AG acceptor sites; such genes generate multiple circular RNAs. For example, DIS3L2, an important factor in mRNA degradation,, regulates Let-7g expression and apoptosis by targeting TGFBR1 and down-regulating the activity of TGF-β pathway. In our study, the DIS3L2 gene was found to yield three circRNA isoforms, ssc-ciR-05526, ssc-ciR-05527, and ssc-ciR-05528, of 349 bp, 523 bp, and 603 bp in length, respectively (Fig. 3B). The circular fractions of ssc-ciR-05526, ssc-ciR-05527, and ssc-ciR-05528 (ratio of the circular isoform relative to all transcripts from the same locus) are 3.8, 23.2, and 17.3%, respectively. Many circRNAs have different acceptors, but share the same donor. For example, host gene RBM25 produces two circRNA isoforms as shown in Fig. 3C. This analysis suggests the existence of alternative and interleaved splicing events, where the same splice sites can participate in multiple forward and backward splicing reactions, either to adjacent exons or to distant skipped exons.

We further assessed the abundance of circRNA isoforms, revealing dissimilar expression profiles among different circRNAs originating from the same host gene. Expression profiles of this type have also been observed by other groups., For example, among the three circRNA isoforms of DIS3L2, ssc-ciR-05527 is highly expressed in the lung, testis, and skeletal muscle at D240, whereas ssc-ciR-05526 and ssc-ciR-05528 were only detected in the testis (Fig. 3D). ssc-ciR-03122, the predominant circRNA from the RBM25 gene, is abundantly expressed in the kidney, while its other four circRNA isoforms are weakly expressed or not detected in most tissues (Fig. 3E). Among circRNA isoforms generated from PTGES3, ssc-ciR-01789 is widely expressed in all tissues and abundantly expressed in the testis and skeletal muscle. Expression of ssc-ciR-01789 in skeletal muscle was age-dependent. However, ssc-ciR-01788 from PTGES3 is weakly expressed only in the kidney and testis (Fig. 3F). Another study revealed that some host genes produce multiple circRNAs with dissimilar expression profiles. Collectively, these results suggest that back-splicing by alternative spliceosomes generates circRNA isoforms and contributes to the diversity and functional complexity of circRNAs.

3.5. S. scrofa circRNAs potentially act as miRNA sponges

A previous study indicated that circRNAs function as miRNA sponges to indirectly regulate mRNA expression., To assess whether pig circRNAs regulate gene transcription by binding to miRNAs, we predicted miRNA binding sites of circRNAs using computational methods. We found 4,436 (74.8%), 4,432 (74.7%), and 2,279 (38.4%) putative miRNA-binding sites on 5,394 circRNAs using the micrnada, RNAhybrid, and microtar programs, respectively (Supplementary Table S4, Fig. 4A). Moreover, we observed a significant positive relationship between the number of miRNA binding sites and circRNA length (Supplementary Fig. S2). To reduce the number of predicted false-positive miRNA targets, only miRNAs considered as targets of circRNAs identified using all three programs were further assessed in our study. We identified 1,842 circRNAs with miRNA binding sites that were predicted by all three prediction methods (Fig. 4B); while most circRNAs were predicted to have no potential binding sites (Supplementary Table S4). The miRNA-binding sites were not equally distributed in the circRNAs, only 441 circRNAs out of 1842 (22.20%) had more than 20 potential biding sites (Fig. 4C).

Figure 4

Analysis for circRNAs with the properties of miRNA sponges. (A) Prediction results of miRNA binding sites for circRNA by mirnada, RNAhybrid, and microtar program, respectively. (B) Overlapping of miRNA binding sites for circRNAs among three methods. (C) Number of miRNA binding sites for circRNAs. (D) A view of interaction between miRNA and circRNAs.

Subsequently, we analysed interactions between differentially expressed circRNAs and predicted target miRNAs. The resulting circRNA–miRNA association network provided nodes and connections between circRNAs and their target miRNAs. According to our data, no miRNA was regulated by only one circRNA. For example, miRNA-378 is involved in a regulatory pathway for circRNA ciR-90. In total, 141 miRNAs and 81 circRNA host genes were involved in this network (Fig. 4D). These results suggest that circRNAs could regulate gene expression by functioning as miRNA sponges. We found that each circRNA can contain more than one miRNA binding site and can bind to more than one miRNA. Although these in silico results should be investigated in vivo, these results illuminate the manner in which circRNAs regulate mRNA translation.

3.6. Dominant circRNAs are expressed in a tissue-specific manner and highly enriched in the testis

To act as miRNA sponges or perform other non-catalytic cellular functions, circRNAs must be expressed at a consequential level within a cell or region of tissue. Here, in order to gain insight into the circRNA landscape across different organs and infer the abundance of each circRNA, we multiplied the circular fraction of each circRNA by the density of RNA-seq reads arising from the cognate gene locus (measured in fragment per kilobase of transcript per million fragments sequenced, or FPKM) (Supplementary Table S5). Most circRNAs are expressed with relatively low abundance., In this study, the average abundance of each S. scrofa circRNA was 0.45 FPKM. When considering the 395 circRNAs with inferred FPKM ≥ 1.0, only 43 had FPKM ≥ 10.0, while 4 had FPKM ≥ 100. To test whether circRNAs are expressed in a tissue-specific manner in pigs, a clustered heatmap was generated for nine organs (Fig. 5A and Supplementary Fig. S3). These results suggest that pig circRNAs have highly tissue-specific expression patterns (Fig. 5A). The largest cluster of circRNAs was in the testis and appeared to be testis-specific (Fig. 5A). Interestingly, circRNA ssc-ciR-03062, which is generated from MYH7 and MYH6, which are associated with familial hypertrophic cardiomyopathy,, was the most abundantly expressed circRNA in the heart and was highly expressed in skeletal muscle.

Figure 5

The tissue-specific expression of S. scrofa circRNAs in pigs. (A) Hierarchical clustering of the circular fraction of 5,399 circRNAs expressed in nine tissues with ≥2 reads. (B) Numbers of circRNAs identified each tissues. (C) Number of S. scrofa circRNAs in tissues analysed. (D) Percentage of tissues specific, -enrichment, and -enhanced circRNAs in each tissues. (E) KEGG term enrichments of circRNAs reveal many genes sets that were specific in tissues.

We found that 339–3,169 circRNAs were detected across most tissues types, and the testis had more than nine times as many circRNAs as the spleen (Fig. 5B). Some circRNAs are ubiquitously expressed, but most are found in only a few tissue types (Fig. 5C). We found that only 37 circRNAs are expressed in all nine organs analysed in this study. This variation was not explained by differences in sequencing depth. It is notable that 53.86% of the identified circRNAs (3,196/5,934 = 53.86%) were detected in the testis; however, 70.86% (2,265/3,196 = 70.86%) of the identified circRNAs were only observed in one other tissue type (Fig. 5B). Tissue-specific circRNAs are a group of circRNAs whose function and expression are confined to one or several tissues., As shown in Supplementary Table S6, the testis shows the largest number of tissue-specific circRNAs (1,155), followed by the heart (205 circRNAs), muscle (174 circRNAs), and fat (147 circRNAs). The lowest number of tissue-specific circRNAs was found in the liver (53 circRNAs) (Table 1). It is notable that we did not detect any spleen-specific circRNAs. We assessed enrichment of gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) terms among the tissue-specific circRNAs. In the heart, GO terms related to muscle development, construction, cell differentiation, and actin filament-based processes are highly enriched (Supplementary Table S7). Pathway analysis indicates that the set of circRNAs specific to the heart is significantly enriched in genes related to hypertrophic cardiomyopathy (HCM), cardiac muscle contraction, dilated cardiomyopathy, and the Wnt signalling pathway (Fig. 5E).

Table 1

Classification of circRNA expression in different tissues

Tissues	Tissues-specific circRNAs	Tissues-enrichment circRNAs	Tissues-enhanced circRNAs
Spleen	0	49	182
Lung	130	125	568
Kidney	99	106	433
Liver	53	53	232
Fat	147	149	371
Testis	1,155	1,192	1,596
Muscle	174	186	412
Ovarium	94	95	555
Heart	205	223	512

Furthermore, we also analysed tissue-enriched (expression level at least 5-fold higher in a particular tissue as compared to all other tissues) and tissue-enhanced (at least 5-fold higher expression level in a particular tissue as compared to the average level in all tissues) circRNA. We identified 2,167 tissue-enriched circRNAs (Supplementary Table S8) and 3,301 tissue-enhanced circRNAs (Supplementary Table S9) across all nine assessed tissue types. The numbers of tissue-enriched and tissue-enhanced genes are highly variable among the analysed tissue types (Table 1, Fig. 5D). The testis shows the largest number of tissue-enriched circRNAs (1,192), followed by the heart (223) and skeletal muscle (186). These results indicate that circRNAs are expressed with highly tissue-specific patterns in pigs. A previous study also reported that dominant circular RNAs were expressed in a tissue-specific manner during human fetal development.

3.7. Differentially expressed circRNAs in postnatal skeletal muscle

Skeletal muscles allow movement of the body via energy expenditure. Postnatal growth of skeletal muscle is mainly realized through increases in the length and girth of muscle fibres, but not by increases in muscle fibre number. To understand the regulation of circRNAs in skeletal muscle, we performed expression profiling of circRNAs in skeletal muscle at postnatal days 0, 30, and 240 (designated D0, D30, and D240, respectively). This profiling allowed evaluation of dynamic changes in circRNA expression from birth until adulthood and identification of circRNAs associated with muscle growth. We detected 1,489 circRNAs expressed with more than two reads in at least one sample (Supplementary Table S10). Of these 1,489 skeletal muscle circRNAs, 597 circRNAs were detected at D0, whereas 670 were detected at D30 and 673 were detected at D240; 57.2–63.9% of skeletal muscle circRNAs were detected at only one stage (Fig. 6A). The identified skeletal muscle circRNAs were studied further to assess changes in their expression levels during skeletal muscle growth. To obtain an overview of spatial circRNA expression changes in developmental skeletal muscle, we clustered pig circRNA expression data from D0, D30, and D240, revealing differential circRNA expression between growth stages. Interestingly, cluster analysis showed that skeletal muscle at D30 had a circRNA expression profile more similar to that of skeletal muscle at D240 than to that of skeletal muscle at D0 (Fig. 6B).

Figure 6

Temporal expression of circRNAs in skeletal muscle. (A) Venn chart of circRNAs detected in skeletal muscle. (B) Heatmap of circRNA expression abundance at 0 day, 30 days, and 240 days. (C) Top circRNAs expressed in skeletal muscle. (D–E) Differentially expressed circRNA between different growth stages in postnatal skeletal muscle.

Next, we focused on abundant circRNAs in skeletal muscle. Notably, we found that several of the most abundant circRNAs in skeletal muscle originate from protein coding genes with pivotal roles in skeletal muscle growth and hypertrophy (e.g. MYH1, MYH2, MYH3, MYH6, MYH7, MYL12B, IGFBP5) or non-coding RNAs with roles in slow-twitch fibre formation (miRNA-208b encoded by introns of myosin genes) (Fig. 6C).,MYH1, MYH6, and MYH7 belong to the MYH myosin superfamily. Myosin, the primary component of thick filaments, is widely expressed in mammalian skeletal muscle and involved in muscle contraction, phagocytosis, cell motility, and vesicle transport., Recent studies reported that myosins are associated with skeletal muscle development and influence the quality of pork. We analysed the expression patterns of all circRNAs detected in skeletal muscle, revealing 16 expression patterns and significant enrichment in cluster 5 (Supplementary Fig. S3). GO analysis of circRNAs expressed in skeletal muscle did not reveal significant enrichment in GO terms.

To study functional transitions in skeletal muscle development, we assessed two types of differentially expressed circRNAs (DEC) between time points. We focused on circRNAs with fold-change >2, P-value < 0.05, and FDR < 0.1. The expression levels of most circRNAs are not significantly altered during postnatal skeletal muscle growth. Type 1 DECs were defined as those that were differentially expressed between two closed time points. We detected 101 and 83 differently expressed circRNAs in the D0 versus D30 and D30 versus D240 analyses, respectively. Among these DECs, 75 circRNAs were up-regulated at D30 in comparison with D0, whereas 26 circRNAs were up-regulated at D0 in comparison with D30 (Supplementary Table S11, Fig. 6D). In the D30 versus D240 group, 58 circRNAs were up-regulated at D240 in comparison with D30, whereas 25 circRNAs were up-regulated at D30 in comparison with D240 (Supplementary Table S12, Fig. 6E). Under the assumption that the functions of circRNAs are related to the functions of their host genes, we performed GO and KEGG pathway analysis to predict the functions of DECs during postnatal growth of skeletal muscle. The GO analysis indicated that the host genes of the D0 versus D30 DECs are significantly associated with muscle development (Supplementary Fig. S4A). However, for the D30 versus D240 DECs, host genes are associated with homeostatic processes such as ion homeostasis and cation homeostasis (Supplementary Fig. 4B). Pathway analysis indicated that the host genes of the D0 versus D30 DECs are involved in tight junction signalling (Table 2). In the D30 versus D240 group, circRNA host genes are significantly associated with the calcium signalling pathway and aldosterone-regulated sodium reabsorption signalling pathway. The overall induction of circRNAs during skeletal muscle development indicates that circularization is likely important for muscle function, because many circRNAs are up- or down-regulated at different time points of skeletal muscle growth.

Table 2

Pathway analysis of parent genes of circRNAs

Group	Term	P value	Benjamini q value
D0 versus D30	Tight junction	7.67E−05	0.0029
D30 versus D240	Calcium signalling pathway	0.0365	0.67
D30 versus D240	Aldosterone-regulated sodium reabsorption	0.0703	0.66

Note: Significantly over-represented KEGG (Kyoto Encyclopaedia of Genes and Genomes) pathways detected with DAVID for parent genes of circRNAs that are differently expressed with at least a twofold changes in group D0 versus D30 and D30 versus D240, respectively.

Type 2 DECs were defined as circRNAs that were differentially expressed only in the D0 versus D30 or D30 versus D240 analyses, representing circRNAs that were up- or down-regulated only between two consecutive time points during muscle growth. We found 29 circRNAs that were differentially expressed in the D0 versus D30 comparison, but not in the D30 versus D240 comparison (Supplementary Table 13). Notably, several circRNAs generated from MYH1, MYH2, MYH6, MYH7, and TNNC2, as well as miRNA-208b, were significantly differentially expressed only in the D0 versus D30 comparison. These host genes mainly encoded proteins associated with muscle development and fibre twitching. We found 24 circRNAs that were differentially expressed in the comparison of D30 and D240, but not in the comparison of D0 and D30 (Supplementary Table S14). Interestingly, these circRNAs were mainly generated from NDFIP2, RBM6, RYR1, SEC24A, RYR1, RBM24, SAR1B, PLEKHF1, ATP1A2, ATP1A4, DCUN1D1, ATP1A2, PRSS22, and PHKB, which encode proteins associated with glycosaminoglycan metabolism and calcium channels. These circRNAs could signify the onset/termination of growth and/or physiological processes at a particular developmental stage.

3.8. S. scrofa circRNA database

In order to facilitate dissemination of information regarding circRNA annotation, circRNA expression, potential circRNA functions, and circRNA–miRNA–gene interactions, we created a database of S. scrofa circRNAs that was designated pigcirNet (http://lnc.rnanet.org/circ). The web interface of pigcirNet is summarized in Fig. 7. PigcirNet shows circRNA location, circular nucleotide length, strand signal, gene Entrez ID, Ensembl id, gene symbol, classification, and expression level across all organs and stages selected in this study. Information regarding prediction of circRNA targeting by miRNAs and circRNA–miRNA interactions are also available in pigcirNet. In order to investigate potential functions of circRNAs, we retrieve GO information with host-related genes. We also search the PubMed database to retrieve reports related to queried circRNAs. A Google© chart is used to display a diagram of circRNA abundance, whereas sigma.js is used to show a circRNA, mRNA, and miRNA interaction network. circRNAs can be browsed by their chromosomal location, classification, or tissue-specific expression pattern. We also provide a search box, in which users can find related circRNAs by searching circRNA id, gene symbol, gene Entrez ID, or gene Ensembl ID. To the best of our knowledge, our database is the first public circRNA resource for a non-rodent mammal.

Figure 7

Framework of the database construction in pigcirNet. The graph illustrates how the network in pigcirNet was constructed.

4. Discussion

Genome-wide analyses of RNA-seq data have revealed that circRNAs are abundant in animal transcriptomes and identified thousands of circRNAs in humans,, mice, nematodes, and drosophila., Guo et al. detected 7,112 human circRNAs using 39 biological samples of whole-cell non-poly(A)-selected RNA-seq data from the ENCODE project. Gao et al. applied the CIRI algorithm to identify nearly 98,526 circRNAs based on ENCODE RNA-seq data from 15 cell lines. Based on 464 RNA-seq samples across 26 human tissues and 104 disease conditions, Liu and colleagues identified 212,950 human circRNAs (53,687 novel circRNAs). These findings suggest that a large number of circRNAs remain unknown and indicate that additional circRNAs can be identified by improving RNA analysis methods and accumulating deeper sequencing data. Because circRNAs are expressed in a highly spatio-temporal specific manner, it is essential that studies of circRNAs in mammals assess various tissues, conditions, and developmental stages. The pig is an important farm animal that provides protein for humans and an important non-rodent animal model that is widely used in biomedical research. Morten et al., the first group to study S. scrofa circRNAs, identified 4,634 unique circRNAs from 2,195 host genes in five brain tissue types at six time-points during fetal porcine development, and the results suggested that some circRNAs were conserved between humans and mice. To perform genome-wide identification of S. scrofa circRNAs and explore their spatio-temporal expression patterns, we carried out total RNA sequencing across nine organs (heart, liver, spleen, lung, kidney, ovarium, testis, skeletal muscle, and fat) and skeletal muscle at three developmental stages (0, 30, and 240 days after birth). We identified 5,934 S. scrofa circRNAs, which represent a significant addition to the growing catalogue of annotated mammalian circRNAs. The circRNAs CDR1as and Sry were well-known studied. However, we did not detect circRNAs formed by a CDR1as orthologue or Sry in our data. Sry is induced during embryonic development. Because circRNAs are expressed in a tissue- and cell-specific manner, S. scrofa circRNAs orthologous to CDR1as and Sry were probably not expressed in the tissues selected in our study. To validate this conjecture, we further analysed expression of CDR1 and Sry. As expected, we did not detect expression of CDR1 or Sry in any tissue, Neither CDR1 nor Sry were mentioned in Morten’s study of the pig brain, indicating that these circRNAs were not expressed in these samples. We analysed a limited set of tissues and developmental stages in this study. Thus, the true number of S. scrofa circRNAs is almost certainly much larger than that reported in this study. In comparison with humans, knowledge of S. scrofa circRNAs remains limited. Therefore, deeper RNA-sequencing to identify circRNAs across additional organs and developmental stages in pigs is essential.

Although the molecular functions of circRNAs are mostly unclear, some circRNAs affect gene expression by acting as microRNA sponges. The most emblematical examples of circRNAs that act as microRNA sponges are CDR1as/ciRS-7 and Sry, which function in human and mouse brain development., Guo et al. predicted that circRNAs from the human C2H2 zinc finger gene family function as miRNA sponges. Another study showed that Drosophila circular RNAs harbour >1000 well-conserved canonical miRNA seed matches. To assess whether S. scrofa circRNAs affect post-transcriptional gene regulation by binding to miRNAs, we predicted circRNA-originating targets by three computational methods. We found that 1,842 of 5,934 (31.04%) circRNAs had putative miRNA-binding sites. In a study of rice, 17.33% (235/1356) of circRNAs were found to have putative miRNA-binding sites. The identified circRNAs had, on average, 19 miRNA-binding sites. There was a significant positive relationship between the number of miRNA-binding sites and circRNA sequence length (Supplementary Fig. S2). Based on our analysis of circRNAs, we constructed an interaction network of miRNAs and circular RNAs. Certainly, the potential genome-wide interplay between miRNAs and circular RNAs warrants further experimental and computational investigation in the future. Although the aforementioned studies indicated that some circRNAs function as miRNA sponges, 4,092 of 5,934 (68.96%) S. scrofa circRNAs identified in this study have no putative miRNA binding site. Several studies have suggested that most circRNAs do not act as miRNA sponges., These findings indicate that circRNAs have functions other than binding to miRNAs, including regulation of host gene transcription, protein binding, and translation.,

Spatio-temporal gene expression patterns reveal the functions of genes within specific tissues/cells at specific times during development. Several studies have revealed that circRNAs tend to show tissue-/stage-specific expression and enrichment in the brain. We also identified tissue-specific, tissue-enriched, and tissue-enhanced circRNAs, providing important information regarding their functions. Of the identified heart-specific circRNAs, for example, ssc-ciR-03925, ssc-ciR-03394, ssc-ciR-03072, ssc-ciR-05708, ssc-ciR-02839, ssc-ciR-02761, ssc-ciR-04971, ssc-ciR-04972, and ssc-ciR-03852 were specifically expressed in the heart. ssc-ciR-04971 is generated from host gene RYR2, mutations in which are associated with stress-induced polymorphic ventricular tachycardia and arrhythmogenic right ventricular dysplasia. Ssc-ciR-03072 is produced from host gene MYH7, changes in the relative abundance of which are correlated with the contractile velocity of cardiac muscle; mutations in MYH7 are associated with familial hypertrophic cardiomyopathy, myosin storage myopathy, dilated cardiomyopathy, and Laing early onset distal myopathy., These findings suggest that heart-specific circRNAs might be involved in heart development and cardiac disease.

Skeletal muscle mass increases during postnatal animal development via hypertrophy. The role of circRNAs in skeletal muscle development is unknown. Therefore, we analysed dynamic changes in circRNA abundance during postnatal growth of skeletal muscle. We identified 149 circRNAs that were differentially expressed between the selected time points and thus might be associated with postnatal growth of skeletal muscle. No studies have assessed the role of circRNAs in myogenesis. The D0 versus D30 DECs and D30 versus D240 DECs differed significantly with respect to enriched GO terms and signalling pathways. The D0 vs. D30 DECs were markedly enriched in GO terms associated with muscle contraction, muscle organ development, and chromatin modification, while the D30 versus D240 DECs were enriched in GO terms associated with cellular cation homeostasis, ATP hydrolysis-coupled proton transport, energy-coupled proton transport against the electrochemical gradient, and response to hypoxia. KEGG analysis suggested that DECs were mainly enriched in the tight junction pathway for the D0 versus D30 group, while DECs in the D30 versus D240 group were significantly associated with calcium signalling and aldosterone-regulated sodium reabsorption signalling. Tight junction proteins participate in the regulation of cell proliferation, gene expression, and cell differentiation., A previous study showed that proliferation and fusion of satellite cells, leading to an increase in the number of myonuclei, may also regulate muscle growth during the early, but not late, stages of postnatal development. Our study indicated that circRNAs (such as ssc-ciR-02753, ssc-ciR-04353, ssc-ciR-04335, ssc-ciR-04349, ssc-ciR-04348, ssc-ciR-04359, ssc-ciR-03066, ssc-ciR-03069, and ssc-ciR-03065) may regulate muscle growth by affecting cell proliferation and fusion during early postnatal muscle development. Clustering analysis suggested that the circRNA expression profile of skeletal muscle at D30 was more similar to that of skeletal muscle at D240 than to that of skeletal muscle at D0. We speculate that circRNAs such as ssc-ciR-01595, ssc-ciR-01592, ssc-ciR-03395, ssc-ciR-02589, ssc-ciR-02611, and ssc-ciR-02610 contribute to muscle plasticity and contraction. Some circRNAs are evolutionarily conserved in terms of sequence and expression, suggesting that they possess similar cellular functions in diverse species.,, Our genome-wide analysis of the spatio-temporal expression patterns of S. scrofa circRNAs provides a foundation for studies aimed at understanding the molecular functions of mammalian circRNAs.

Thousands of circRNAs have been discovered in plants, animals, and humans. Given the emerging understanding of the biological importance of circRNAs and research efforts to understand them, researchers have constructed several circRNA databases, including Circ2Traits, CircNet, deepBase v2.0, CircInteractome, and SomamiR2.0. At present, these databases include human, mouse, and C. elegans circRNAs, but not S. scrofa circRNAs, because identification of circRNAs in pigs has lagged far behind identification of circRNAs in humans, mice, and C. elegans. Exploring the pig genome can illuminate facets of human culture, organic evolution, biomedical research, and animal breeding. Therefore, we constructed the database of S. scrofa circRNAs based on our data. To the best of our knowledge, this database is the first public resource containing information about the circRNAs of a non-rodent mammal.

Availability

RNA-seq data were deposited in the Gene Expression Omnibus with accession codes was GSE73763. The identified circRNA was stored in http://lnc.rnanet.org/circ.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.