Skip to main content
Download PDF

Genome-wide analysis of circRNA regulation during spleen development of Chinese indigenous breed Meishan pigs

BMC Genomics volume 24, Article number: 477 (2023) Cite this article

Abstract

Background

Numerous circular RNAs (circRNAs) have been recently identified in porcine tissues and cell types. Nevertheless, their significance in porcine spleen development is yet unelucidated. Herein, we reported an extensive overlook of circRNA expression profile during spleen development in Meishan pigs.

Results

Overall, 39,641 circRNAs were identified from 6,914 host genes. Among them, many circRNAs are up- or down-regulated at different time points of pig spleen development. Using WGCNA analysis, we revealed two essential modules for protein-coding genes and circRNAs. Subsequent correlation analysis revealed 67 circRNAs/co-expressed genes that participated in immnue-associated networks. Furthermore, a competing endogenous RNA (ceRNA) network analysis of circRNAs revealed that 12 circRNAs modulated CD226, MBD2, SAMD3, SIT1, SRP14, SYTL3 gene expressions via acting as miRNA sponges. Moreover, the circRNA_21767/miR-202-3p axis regulated SIT1 expression in a ceRNA manner, which is critical for the immune-based regulation of spleen development in Meishan pigs.

Conclusions

Overall, our results demonstrated that the circRNAs were differentially expressed during different stages of porcine spleen development, meanwhile the circRNAs interacted with immune-related genes in a ceRNA-based fashion. Moreover, we presented biomedical researchers with RNAseqTools, a user-friendly and powerful software for the visualization of transcriptome profile data.

Peer Review reports

Introduction

Pigs (Sus scrofa) are an incredible resource for agriculture and food production, as well as biomedical research for the examination of human development, congenital diseases, and pathogen response networks [1]. Till date, systemic research facilitated the pig industry to critically select for enhanced growth, feed conversion efficiency, and disease resistance. In fact, disease resistant pig breeds like Yorkshire and Landrace experience considerably fewer diseases than the local indigenous pigs in China [2, 3]. The Meishan pigs is one of four popular indigenous breeds in China, with a history spanning over 400 years. It is particularly known for its enhanced fertility, strong immune response, and excellent meat quality. Moreover, some investigations reported that the Chinese native-bred Meishan pigs may possess stronger immunity and disease resistance than the commercial or crossbred pigs [4, 5]. Nevertheless, the genetic basis and regulation mechanism of the immune system are still not fully understood and need to be further explored.

Circular RNA (circRNA) is a newly discovered form of non-coding RNA (ncRNA) known for its enclosed configuration, generated via precursor mRNA back-splicing [6]. In contrast to linear RNAs, circRNAs lack both the 5′ cap and 3′ tail, and are covalently closed, thereby making them more stable, RNase R resistant, and with prolonged half-lives [7]. Given these unique properties, circRNAs are great potential as an optimal biomarker. With increasing advancement in high-throughput RNA-seq technology, a growing number of studies identified circRNA presence and essential roles in various organisms. Emerging evidences revealed that circRNAs critically modulate several physiological processes [8, 9], including immune cells and various immune responses [10, 11]. CircRNAs act in a tissue and developmental stage-specific manner, and multiple reports confirmed their functions in the porcine skeletal muscle [12], embryonic muscle [13], Liver [14], longissimus thoracis [15], brain [16], etc. Unfortunately, studies involving the growth and development of porcine immune organs is incredibly rare. Thus, it is challenging to fully elucidate the modulatory properties of circRNAs in immune cells and immune responses. Hence, herein, we attempted to explore the circRNAs-mediated regulation networks that influence porcine immune organs development.

Spleen is a secondary lymphoid organ involved in older erythrocyte elimination, iron recycle, capture and destruction of pathogens, and induction of adaptive immune responses [17, 18]. Being a critical immunologic organ, there is a strong need to examine spleen development in Meishan pigs. Herein, we systematically identified and characterized circRNAs within the porcine spleen from 8 developmental stages, namely, 1d, 7d, 14d, 21d, 28d, 35d, 120d, and 180d after birth, using an Illumina HiSeq platform. Then, using the weighted gene co-expression (WGCNA) and competing endogenous RNA (ceRNA) interaction network analyses, we demonstrated a critical role of circRNAs during porcine spleen development. Our findings provide a beneficial resource for the circRNAs-mediated modulatory functions in Meishan pigs, as well as annotation of the porcine genome. We also contribute to the enhanced understanding of the mammalian spleen development.

Results

Overview of the sequencing information

To explore the presence of circRNAs during spleen development, we assessed circRNAs expression in the spleen tissues of Meishan pigs at various developmental stage. We prepared and sequenced ribo-depleted total RNA-seq libraries, as shown in the flow chart (Fig. 1). Table S2 presents our rudimentary sequencing results. At various developmental stages, we acquired around 111.37 million raw reads and 108.68 million clean reads. Additionally, the error rate of all 24 sequencing data sets was below 0.03%, with quality base scores of Q30 > 95%. Together, these data validated the authenticity and precision of our sequencing data.

Fig. 1
figure 1

A schematic workflow of circRNAs-mediated regulation of spleen development in the Chinese indigenous breed Meishan pigs. Spleen samples were acquired from 8 distinct developmental stages (1d, 7d, 14d, 21d, 28d, 35d, 120d, and 180d after birth), prior to circRNA sequencing (n = 3 samples per time point), and subsequent downstream analyses

Full size image

Identification and characterization of circRNAs expression in Meishan pigs

To characterize the circRNAs in the Meishan pig spleen tissues, we evaluated their sequences. In all, 39,641 candidate circRNAs were identified from 6,914 host genes, using ≥ 1 unique back-spliced read (Fig. 2A). This indicated that most genes potentially produced numerous circRNAs. Their genomic loci showed wide distribution across all chromosomes (Fig. 2B), most (38,507, 97.14%) circRNAs were over 200 nt in length (Fig. 2C), and the GC content of most circRNAs were between 35 and 50% (Fig. 2D). Furthermore, the circRNAs belonged to 3 categories, namely, exonic circRNAs, intronic circRNAs and intergenic. In case of sense circRNAs, most circRNAs (37,101, 97.61%) belonged to the exonic circRNAs category, then the intronic circRNAs category (461, 1.21%), and a small proportion belonged to the intergenic_downstream (225, 0.59%) or intergenic_upstream (223, 0.59%) categories (Fig. 2E). Using principal component analysis (PCA), we revealed that the sample repeatabilities of individual groups were satisfactory (Fig. 2F), confirming that these data can be employed for subsequent DE circRNAs analyses.

Fig. 2
figure 2

Identification and characterization of circRNAs within Meishan pig spleen. (A) circRNAs identification according to the GT-AG backsplicing signal. (B) Length allocation of circRNAs. (C) Chromosomal allocation of circRNAs. (D) GC content frequency allocation of circRNAs. (E) circRNAs classification according to the Sus scrofa genome. (F) Principal component analysis of 24 porcine spleen samples

Full size image

DE circRNAs during spleen development of Meishan pigs

Emerging evidences support the presence of DE circRNAs during porcine skeletal muscle development [19], however, only limited number of studies explored their temporal expression during spleen development. Herein, we examined the splenic circRNAs expression profiles during different stages of development, namely 1, 7, 14, 21, 28, 35, 120, and 180 days after birth. With 1d (just after birth) as the control, we observed 308, 224, 262, 397, 397, 564, and 717 DE circRNAs during the 7d vs. 1d, 14d vs. 1d, 21d vs. 1d, 28d vs. 1d, 35d vs. 1d, 120d vs. 1d, 180d vs. 1d groups, respectively (Fig. 3). Taken together, these data suggested that the circular nature of circRNAs was potentially critical for its immunomodulatory activity, as multiple circRNAs were either highly expressed or scarcely expressed during distinct stages of porcine spleen development.

Fig. 3
figure 3

Volcano plot analysis of differentially expressed (DE) circRNAs during different stages of Meishan pig spleen development. Volcano plots depicting –log10 (pvalue) versus log2foldchange in circRNAs expression profile in RPM over different time periods. Red dots indicate highly expressed circRNAs; Green dots indicate scarcely expressed circRNAs.

Full size image

Generation of the circRNA-mRNA co-expression network using WGCNA

WGCNA is often used to generate gene co-expression modules from gene expression profiles [20]. Herein, we elucidated the circRNAs expression profile by generating a co-expression module containing 3000 circRNAs using WGCNA. With a soft threshold power β of 10 (Fig. 4A), our generated network was categorized into 7 modules, marked by different colors (black, brown, grey, pink, red, turquoise, and yellow) (Fig. 4B, Table S3). Figure 4 C illustrates the association between co-expression modules and various spleen developmental stages. We observed a strong direct correlation between the red modules and developmental stages (correlation coefficient = 0.88, p = 2e-08). The heatmap shows the expression pattern of red module, where the expression of genes in the red module tends to increase with time (Fig. 4D). Hence, the red module for circRNA expression profile was selected for subsequent analysis.

Fig. 4
figure 4

WGCNA analysis of the circRNA dataset identified by the co-expression module of Meishan pig spleen development. (A) Scale independence and average connectivity analysis for varying soft cut-off powers. (B) Clustering dendrograms of porcine spleen samples at various developmental stages. Differing colors represent WGCNA-generated co-expression modules. Branch modules are color-coded based on their interconnections with other circRNAs. Overall, 7 modules of various colors are presented in the horizontal bar using 0.25 threshold merging. (C) The module–trait association. Individual rows represent a modular eigengene, and individual columns represent a trait. Each cell incorporates associated correlation and p-value. Red refers to direct associations, and blue denotes inverse associations. (D) Heatmap depicts the circRNA expression profile of the red module

Full size image

We conducted further analysis on a WGCNA module consisting of 15,188 mRNAs, based on a prior transcriptome mRNA profiling of the Meishan pig spleen at distinct developmental stages (NCBI SRA repository under the BioProject IDs: PRJNA875155, PRJNA869328). As depicted in Fig. S1 and Table S4, based on the soft threshold power β = 9 (Fig. S1A), we generated a gene hierarchy clustering tree, and screened 12 gene modules (brown, darkmagenta, darkolivegreen, darkturquoise, sienna3, lightgreen, greenyellow, skyblue, paleturquoise, lightcyan, royalblue, and grey) that were intricately linked to porcine spleen development (Fig. S1B). Additionally, our assessment of the module-trait association revealed that distinct spleen developmental stages were strongly and positively associated with the darkturquoise module (correlation coefficient = 0.74, p = 4e-05) (Fig. S1C). Moreover, our heatmap revealed the transcript expression profile of darkturquoise (Fig. S1D). Therefore, we employed the darkturquoise module of mRNAs expression profile for further analyses. Collectively, our circRNA-mRNA co-expression modules revealed distinct expression patterns, and preliminary screening revealed that the circRNAs in the red module and mRNAs in the darkturquoise module displayed a strong direct correlation with spleen developmental stage, thereby confirming a significant role of circRNAs in pig spleen development.

Based on the above results, we identified 36 circRNAs in the red module, and 508 mRNAs in the darkturquoise module. We performed correlation analysis on the circRNAs/mRNAs co-expression (Fig. 5A). Our results revealed 13 circRNAs and 54 mRNAs that were strongly associated in their expression profiles, carrying a Pearson correlation coefficient > 0.85 (Fig. 5B, Table S5). We next conducted GO terms and KEGG network annotation analyses of the circRNA-associated genes. Using GO analysis, we revealed that the associated genes participated in immune regulation (Fig. 5C, Table S6). Using KEGG network analysis, we revealed that the T cell receptor signaling pathway, antigen processing and presentation, primary immunodeficiency, intestinal immune network for IgA production were the most immune-associated networks (Fig. 5D, Table S7).

Fig. 5
figure 5

Correlation analysis and functional annotation of circRNA and mRNA expression profiles. (A) Correlation analysis of the circRNAs/mRNAs co-expression profile. (B) Heatmap of circRNA and mRNA gene expression profiles. (C) GO functional analysis of mRNAs associated with the DE circRNAs. (D) KEGG network analysis of mRNAs associated with the DE circRNAs.

Full size image

Generation of a possible circRNA-miRNA-mRNA modulatory axis

To further elucidate the influence of circRNAs on mRNAs, in presence of miRNAs, during Meishan pig spleen development, we generated a ceRNA network using the aforementioned data. Table S8 illustrates the miRNA target gene prediction via the miRanda and Pita software. In all, we utilized 12 circRNAs, 24 miRNAs, and 6 mRNAs to generate the ceRNA network (Fig. 6A). Subsequently, we employed a heatmap to illustrate the expression profiles of CD226, MBD2, SAMD3, SIT1, SRP14, SYTL3, and their co-expressed circRNAs, with a rising trend over increasing developmental stage (Fig. 6B). We also screened for the major immune-related spleen networks, and discovered the SIT1 gene to be a marked contributor of immune regulation (Fig. 6C). Additional correlation analysis revealed that all 9 circRNA expressions were strongly and directly associated with SIT1 expression (Fig. 7). Collectively, these results suggested that circRNAs were DE across all 8 analyzed stages of porcine spleen development. Therefore, they can potentially modulate spleen development via regulation of host gene expression.

Fig. 6
figure 6

Generation of a potential circRNA-miRNA-mRNA modulatory axis. (A) Sankey plot decpiting the ceRNA axis estimated by the targeted associaiton and co-expression of circRNA-miRNA and miRNA-mRNA. Individual rectangles represent a single gene, and the connection degrees of individual genes are visualized according to the rectangular size. (B) Heat-map depicting the pattern of gene and circRNAs expression within the co-expression network. (C) Functional annotation of major genes within the co-expression network

Full size image
Fig. 7
figure 7

Scatter plot of circRNAs and mRNAs correlation analysis during Meishan pig spleen development. Individual genes and circRNAs depicting various levels of correlation. R refers to the correlation coefficient; P < 0.01 denotes significant correlation

Full size image

Experimental verification of the circRNA-miRNA-mRNA regulatory network

Emerging evidences suggested that a series of sub-graphs containing 3 nodes modulate several biological functions [21]. Using circRNA network analysis, we identified a regulatory network with circRNA_21767, miR-202-3p, and SIT1 (Fig. 8A). To confirm this network, we further evaluated the circRNA_21767 and SIT1 contents using RT-qPCR at different points of spleen development. We demonstrated a close association between circRNA_21767 and SIT1 expressions, suggesting that circRNA_21767 potentially modulated SIT1 expression (Fig. 8B C). To further validate the association between miR-202-3p, circRNA_21767 and SIT1, we identified the potential miR-202-3p docking site within circRNA_21767 and 3’UTR of SIT1 (Fig. 8D-F), using dual luciferase reporter assay. We demonstrated that miR-202-3p strongly suppressed circRNA_21767 luciferase activity. Moreover, miR-202-3p considerably diminished SIT1 luciferase activity (Fig. 8G-H), relative to NCs (miR-NC and miR-202-3p inhibitor). Based on these evidences, circRNA_21767 potentially served as an miR-202-3p sponge to regulate SIT1 expression. Taken together, these results revealed that three nodes circRNA mediated regulatory circuitry might play an essential role during spleen development in Meishan pigs.

Fig. 8
figure 8

Experimental validation of circRNA-miRNA-mRNA regulatory circuitry. (A) Illustration of a three nodes circuitry including miR-202-3p, circRNA_21767 and SIT1. RT-qPCR shows the relative expression level of circRNA_21767 (B), miR-202-3p (C) and SIT1 (D) during pig spleen development. (E) Predicted miR-202-3p binding site on circRNA_21767. The design of luciferase reporter. WT, the WT sequence of circRNA_21767 contains miR-202-3p binding site; Mut, the sequence of circRNA_21767 with mutation in miR-202-3p binding site. (F) PK15 cells were co-transfected with wild-type (WT) or mutant (MUT) luciferase reporters of circRNA_21767 with miR-202-3p mimics or negative control (NC) mimics. The relative levels of firefly luminescence normalized to Renilla luminescence are plotted. (G-H) Predicted miR-30a-3p binding site on SIT1-3’UTR. PK15 cells were co-transfected with wild-type (WT) or mutant (MUT) luciferase reporters of SIT1 with miR-202-3p mimics, NC mimics. All data are presented as the mean ± SD, **P < 0.01, nsP > 0.05

Full size image

RNAseqTools construction for the visualization of transcriptome data

To assist biomedical researchers who lack programming expertise in analyzing and visualizing transcriptome sequencing data, we developed a user-friendly software named RNAseqTools (https://github.com/ChaoXu1997/RNAseqTool) using the Shiny framework (Fig. S2A). This interactive toolkit provides diverse analytical functions and visualization modes, such as, PCA, multi-group DE analysis, gene functional annotation, volcano plot, GSEA, gene trend analysis, and WGCNA (Fig. S2B). Its interactive visual interface will allow researchers to effortlessly elucidate data characteristics and distribution for basic research purposes.

Discussion

Genome-wide analyses of RNA-seq data have revealed the abundance of circRNAs in animal transcriptomes and identified thousands of circRNAs in various species including humans, mice, livestock, and poultry, et al. [22]. CircRNAs have been identified in large numbers by deep RNA sequencing and are believed to be widely expressed in eukaryotic cells [23]. These findings suggest that there is still a considerable number of circRNAs that have not been discovered. Recent advancements in transcription profiling technology facilitates our growing understanding of non-coding RNA (ncRNA) and its immunomodulatory function. Investigations involving various animal species reported that circRNAs strongly modulate body’s normal immunity and immune imbalance, which is intricately linked to the development and progression of multiple immune diseases [24]. Moreover, circRNAs is also reported to function in a tissue- and pathogenesis-specific manner in many diseases [25]. Given its highly regulated profile of T and B lymphocytes, the spleen is the most essential organ for combating bacterial and fungul infections [18, 26]. Recent reports revealed the circRNAs characteristics across many tissues in pigs [19], as well as the DE gene and long ncRNA expressions during postnatal spleen development [27]. However, the identities and expression profiles of S. scrofa circRNAs remain largely uncharacterized in spleen tissue.

Spatio-temporal gene expression patterns elucidate the functions of genes within specific tissues or cells at particular times throughout the process of development. Existing research demonstrates the spatio-temporal regulation of circRNA expression in Drosophila brains [28], pig brains [16], and mouse pre-implantation embryos [29]. Additionally, circRNAs exhibit tissue-specific expression in the human heart and lung during fetal development [30]. Due to their precise spatio-temporal expression, it is essential that studies of circRNAs in mammals assess various tissues and developmental stages. Comprehensive examination of the spatio-temporal distribution of circRNAs is crucial for understanding their biological functions. The local pig breeds in China have underdone selective breeding for many generations, and have retained genetic diversities in existing pig populations. Hence, the domestic pig is an excellent model for spleen developmental study in mammals [31]. Herein, we performed genome-wide analysis of circRNA in spleen tissue at multiple developmental stages (from birth to adult) from Meishan pigs (a native Chinese pig), and our screening uncovered 39,641 circRNAs in the developing spleen of Chinese indigenous Meishan pigs. Based on our preliminary literature screening, this report is the first to extensively explore the circRNAs expression profile during the various developmental stages of pig spleen using RNA-seq data. Therefore, the conclusions of this study will be highly beneficial in the exploration and elucidation of circRNAs functions during spleen development. Our findings will also provide a theoretical basis for future investigations on human spleen tissue development and dysfunction.

CircRNAs are known to modulate numerous physiological processes using distinct signaling networks [32]. Although the molecular functions of circRNAs are mostly unclear, some circRNAs affect gene expression by acting as microRNA sponges. For example, CDR1as functions as an efficient miRNA sponge, binding miR-7 and allowing mRNAs to escape degradation following miRNA binding [33]. Besides, circRNAs serve as miRNAs “sponges” to control pig muscle development [13, 15, 34, 35]. However, its role in ceRNA regulation during spleen development remains unclear. Herein, we demonstrated that CD226, MBD2, SAMD3, SIT1, SRP14, and SYTL3 were the strong modulators of spleen development, and these genes were further regulated by bundles of circRNAs. Subsequently, using KEGG analyses, we revealed that the cell adhesion molecules (CAMs) and immunomodulation of several genes were significantly enriched. Cell adhesion molecules, which express on the cell surface, facilitate activities of several biological processes, including, immune response, inflammation, and embryogenesis [36]. CD226 is exclusively expressed on the surface of Th1 cells, and they modulate Th1 differentiation [37]. Herein, we revealed that SIT1 was strongly modulated by the circRNA_21767/miR-202-3p axis. Based on a prior investigation, signaling threshold regulating transmembrane adaptor 1 (SIT1), encoded by the SIT1 gene, is a disulfide-linked homodimeric glycoprotein that interacts with the lymphocyte-specific transmembrane adaptor protein family to modulate the immunologic process [38]. It is implied that these genes can act as checkpoints to activate the immune response in spleen, and inhibitors of these genes (such as miRNA or circRNA) may be used for new non-antibody immune checkpoint inhibitors. The findings from this report provide novel evidence regarding circRNAs involvement in pig spleen immunity. Even though our software-simulated conclusions require further validation in vivo, our results confirmed that circRNAs indeed serve as miRNA sponges to modulate immune-associated gene expression within the spleen. Meanwhile, the non-coding RNA results from transcriptome analysis can serve as beginning experimental targets future in vivo verifications in animal models.

Conclusion

In conclusion, this study presented the first reported catalog of circRNAs expression profile in the spleen tissues from Meishan pigs. We characterized numerous DE circRNAs that were relevant to the porcine spleen development. We also generated a comprehensive expression profile of various RNAs, particularly mRNAs and circRNAs, which, together, modulate a myriad of physiological processes within the Meishan pig spleen. Moreover, we identified a ceRNA network involving circRNA_21767/miR-202-3p/SIT1, which may modulate immune activities. Our findings provide candidate ncRNA-based targets for future in vivo animal study on the molecular regulation of immune function within the porcine spleen. Finally, we designed a user-friendly software named RNAseqTools that is available for use for the visualization of transcriptomic data using the Shiny app (https://github.com/ChaoXu1997/RNAseqTool).

Materials and methods

Animal treatment and sample collection

Overall, 24 healthy Meishan pigs were acquired from the Kunshan Conservation Ltd. (Suzhou City, Jiangsu Province, China) (Permission No. JS-C-05). Our examination involved 8 developmental stages, namely, 1d, 7d, 14d, 21d, 28d, 35d, 120d, and 180d after birth. All piglets were provided with the same standard diet, with housing at an environmentally regulated facility. Three piglets were sacrificed per treatment via an intravenous administration of pentobarbital sodium to minimize suffering. The spleen was immediately extracted and frozen in liquid nitrogen. Our animal care and treatment protocols were approved by the Institutional Animal Care and Use Committee (IACUC) of Yangzhou University (No. SYXK(Su)2021-0026).

RNA extraction, quality assessment and library construction

Splenic total RNA extraction employed TRIZOL (Invitrogen, Carlsbad, CA, United States) following kit directions. RNA quality assessment employed an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States). Ribosomal RNA elimination from extracted RNA utilized the Ribozero™ rRNA Removal Kit (Epicenter, United States). Subsequently, linear RNA was eliminated via an RNAse R kit (Epicenter, United States). Lastly, sequencing libraries were generated with rRNA-free and linear RNA-free RNA using the NEBNext® Ultra™ Directional RNA Library Prep Kit (NEB, Ipswich, MA, United States).

CircRNA sequencing and screening

Sequencing of generated libraries was done via the Illumina HiSeq 2500 platform (Oebiotech Corporation, Shanghai, China). TopHat2 (v2.1.1) mapped the clear data against the porcine reference genome (Sscrofa11.1). We uploaded the resulting RNA-seq information in the Gene Expression Omnibus (accession codes GSE228936). To construct the SAM file, we employed BWA [39] for alignment of the sequencing reads of individual samples to the porcine genome. Lastly, we identified circRNA using CIRI2 tools (version 1.2) [40]. The circRNA sequence estimation employed junction reads and GT-AG cleavage signals. The identified circRNAs were further assessed according to their type, chromosome distribution, and length distribution using find_circ (version 1) [41] and the annotation information of the Sscrofa11.1 genome.

Analysis of the differentially expressed (DE) circRNAs

DE circRNAs were identified using the DESeq2 package [42], based on the following criteria: absolute fold change ≥ 2 and adjusted p-value < 0.05.

Generation of WGCNA

We generated a co-expression network across all examined developmental stages using the R package WGCNA [43]. This allowed us to further analyze protein-coding genes and circRNAs sets that were critical at each stage of development.

Establishment of the circRNA–miRNA–mRNA ceRNA axis

We retrieved the miRNA sequences from the miRbase (https://mirbase.org/, accessed on 20 December 2022) website, and the mRNA 3’UTR sequences from BioMart (https://asia.ensembl.org/index.html, downloaded on 25 December 2022) [44]. To explore the associations and functions of various ncRNAs and mRNAs, we next generated an ncRNA–mRNA modulatory axis. MiRanda (version:3.3a) [45] predicted the circRNA–miRNA–mRNA pairs (score threshold ≥ 150 and energy threshold ≤ − 20). The circRNA–miRNA–mRNA co-expression axis was then generated with the Cytoscape software (v3.2.1) [46] to examine the role of major circRNAs. Finally, we visualized the network using the ggalluvial R package (Version: 0.9.1).

Functional enrichment analysis (FEA)

We performed the GO term and KEGG pathway [47,48,49] enrichment analyses of circRNA-associated genes using the R package clusterProfiler [50], using the following criteria: corrected P-value < 0.05 was the significance threshold.

Quantitative RT-PCR (qPCR) validation

qPCR employed an Applied Biosystems qPCR apparatus and SYBR green master mix (Vazyme, Nanjing, China), as per kit protocols. The reaction mixture was composed of 5 µl SYBR Green Mixture (Vazyme, Nanjing, China), 1 µl cDNA template, 0.2 µl primer (each), and 3.6 µl deionised water, and the conditions were set as follows: 95 °C for 5 min, 40 cycles of 95 °C for 10 s, and 60 °C for 30 s. All employed primer sequences were summarized in Table S1. GAPDH (for mRNA and circRNA) and U6 (for miRNA) served as controls for normalization purpose, and all qRT-PCR reactions were performed three times, and the average adjusted value was used for analysis. Relative gene expressions were computed using the 2−ΔΔCt formula. All data analysis employed the unpaired 2-tailed Student’s t-tests; *p < 0.05, **p < 0.01, and data are provided as mean ± SD (standard error of mean).

Dual-luciferase reporter assays

We inserted wild-type (WT) or mutant-type (Mut) SIT1-3’UTR or circRNA_21767 downstream of the pmirGLO Dual-Luciferase vector. Next, we plated cells into a 24-well plates, and co-transfected the cells with the WT (SIT1-WT or circRNA_21767-WT), or Mut (SIT1-Mut, circRNA_21767-Mut), and miR-202-3p-mimics (NC-mimics). Finally, we detected both firefly and renilla bioluminescence using the DLR Gene Assay Kit (Beyotime, Shanghai, China). The firefly luciferase activity was then adjusted with the control renilla luciferase activity for a measurement of the constructed reporter luciferase activity.

SIT1-WT:

TCCTCCCAACCCCCAAACTCCCAGGTTTTCAGTCCTCCTTCCGGAGTTTAATCAGATGCTCCCCACTCCGGCTGCCTCATGGTCTGTGCTGACGCCTCAGTGTCTCCTCAGCCACAGGAAGTAGGCAGTGGGGGAGGGGGTTAGAGCCTGAGAGGATATGTATGGGTATC.

SIT1-Mut:

TCCTCCCAACCCCCAAACTCCCAGGTTTTCAGTCCTCCTTCCGGAGTTTAATCAGATGCTCCCCACTCCGGCTGCCTCATGGCTCACATCGGTATTCTAGTGTCTCCTCAGCCACAGGAAGTAGGCAGTGGGGGAGGGGGTTAGAGCCTGAGAGGATATGTATGGGTATC.

circRNA_21767-WT:

TAGTCTTGTATTGATGGTTTTGCACTATTTACAGAccctACCTGAAcctatccttccatccatccaaaAAATGTACCCAGAATCTTTTAGTCCTGCCCTACAGCTGCACCTCGTACATCAAGCTCCATATAATGTTCCTCCTTACCTCTCAAAGAACGAATCAAATCTTGGGGACCTCTTAT.

circRNA_21767-Mut:

TAGTCTTGTATTGATGGTTTTGCACTATTTACAGAccctACCTGAAcctatccttccatccatccaaaAAATGTATTTGGAATCTTTTAGTCCCATTCTACAGCCATGTTCTGTACATCAAGCTCCATATAATGTTCCTCCTTACCTCTCAAAGAACGAATCAAATCTTGGGGACCTCTTAT.

Data Availability

The circRNA sequencing datasets during the current study are available in the Gene Expression Omnibus repository (accession codes GSE228936).

References

  1. Lunney JK, Van Goor A, Walker KE, Hailstock T, Franklin J, Dai C. Importance of the pig as a human biomedical model. Sci Transl Med. 2021;13:eabd5758.

    CAS  PubMed  Google Scholar 

  2. Clapperton M, Bishop SC, Glass EJ. Innate immune traits differ between Meishan and large White pigs. Vet Immunol Immunopathol. 2005;104(3–4):131–44.

    CAS  PubMed  Google Scholar 

  3. Liang W, Li Z, Wang P, Fan P, Zhang Y, Zhang Q, Wang Y, Xu X, Liu B. Differences of immune responses between Tongcheng (chinese local breed) and large White pigs after artificial infection with highly pathogenic porcine reproductive and respiratory syndrome virus. Virus Res. 2016;215:84–93.

    CAS  PubMed  Google Scholar 

  4. Duchet-Suchaux MF, Bertin AM, Menanteau PS. Susceptibility of chinese Meishan and European large white pigs to enterotoxigenic Escherichia coli strains bearing colonization factor K88, 987P, K99, or F41. Am J Vet Res. 1991;52(1):40–4.

    CAS  PubMed  Google Scholar 

  5. Chen J, Qi S, Guo R, Yu B, Chen D. Different messenger RNA expression for the antimicrobial peptides beta-defensins between Meishan and crossbred pigs. Mol Biol Rep. 2010;37(3):1633–9.

    CAS  PubMed  Google Scholar 

  6. Li X, Yang L, Chen LL. The Biogenesis, Functions, and Challenges of Circular RNAs. Mol Cell. 2018;71(3):428–42.

    CAS  PubMed  Google Scholar 

  7. Chen LL, Yang L. Regulation of circRNA biogenesis. RNA Biol. 2015;12(4):381–8.

    PubMed  PubMed Central  Google Scholar 

  8. Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15(7):409.

    PubMed  PubMed Central  Google Scholar 

  9. Jeck WR, Sharpless NE. Detecting and characterizing circular RNAs. Nat Biotechnol. 2014;32(5):453–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Zhang Y, Zhang Y, Li X, Zhang M, Lv K. Microarray analysis of circular RNA expression patterns in polarized macrophages. Int J Mol Med. 2017;39(2):373–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Huang ZK, Yao FY, Xu JQ, Deng Z, Su RG, Peng YP, Luo Q, Li JM. Microarray expression Profile of Circular RNAs in Peripheral Blood mononuclear cells from active tuberculosis patients. Cell Physiol Biochem. 2018;45(3):1230–40.

    CAS  PubMed  Google Scholar 

  12. Ma L, Chen W, Li S, Qin M, Zeng Y. Identification and functional prediction of circular RNAs related to growth traits and skeletal muscle development in Duroc pigs. Front Genet. 2022;13:858763.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Hong L, Gu T, He Y, Zhou C, Hu Q, Wang X, Zheng E, Huang S, Xu Z, Yang J, et al. Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development. Front Cell Dev Biol. 2019;7:289.

    PubMed  PubMed Central  Google Scholar 

  14. Chen W, Ma H, Li B, Yang F, Xiao Y, Gong Y, Li Z, Li T, Zeng Q, Xu K, et al. Spatiotemporal regulation of circular RNA expression during Liver Development of Chinese Indigenous Ningxiang Pigs. Genes (Basel). 2022;13(5):746.

    CAS  PubMed  Google Scholar 

  15. Zhuang X, Lin Z, Xie F, Luo J, Chen T, Xi Q, Zhang Y, Sun J. Identification of circRNA-associated ceRNA networks using longissimus thoracis of pigs of different breeds and growth stages. BMC Genomics. 2022;23(1):294.

    CAS  PubMed  PubMed Central  Google Scholar 

  16. Venø MT, Hansen TB, Venø ST, Clausen BH, Grebing M, Finsen B, Holm IE, Kjems J. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol. 2015;16:245.

    PubMed  PubMed Central  Google Scholar 

  17. Bronte V, Pittet MJ. The spleen in local and systemic regulation of immunity. Immunity. 2013;39(5):806–18.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Lewis SM, Williams A, Eisenbarth SC. Structure and function of the immune system in the spleen. Sci Immunol. 2019;4(33):eaau6085.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Liang G, Yang Y, Niu G, Tang Z, Li K. Genome-wide profiling of Sus scrofa circular RNAs across nine organs and three developmental stages. DNA Res. 2017;24(5):523–35.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:Article17.

    PubMed  Google Scholar 

  21. Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, Alon U. Network motifs: simple building blocks of complex networks. Science. 2002;298(5594):824–7.

    CAS  PubMed  Google Scholar 

  22. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495:333–8.

    CAS  PubMed  Google Scholar 

  23. Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L. Complementary sequence-mediated exon circularization. Cell. 2014;159:134–47.

    CAS  PubMed  Google Scholar 

  24. Geng H, Tan XD. Functional diversity of long non-coding RNAs in immune regulation. Genes Dis. 2016;3(1):72–81.

    PubMed  PubMed Central  Google Scholar 

  25. Ng WL, Marinov GK, Liau ES, Lam YL, Lim YY, Ea CK. Inducible RasGEF1B circular RNA is a positive regulator of ICAM-1 in the TLR4/LPS pathway. RNA Biol. 2016;13(9):861–71.

    PubMed  PubMed Central  Google Scholar 

  26. Mebius RE, Kraal G. Structure and function of the spleen. Nat Rev Immunol. 2005;5(8):606–16.

    CAS  PubMed  Google Scholar 

  27. Che T, Li D, Jin L, Fu Y, Liu Y, Liu P, Wang Y, Tang Q, Ma J, Wang X, et al. Long non-coding RNAs and mRNAs profiling during spleen development in pig. PLoS ONE. 2018;13(3):e0193552.

    PubMed  PubMed Central  Google Scholar 

  28. Westholm JO, Miura P, Olson S, Shenker S, Joseph B, Sanfilippo P, Celniker SE, Graveley BR, Lai EC. Genome-wide analysis of Drosophila circular RNAs reveals their structural and sequence properties and age-dependent neural accumulation. Cell Rep. 2014;9(5):1966–80.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Fan X, Zhang X, Wu X, Guo H, Hu Y, Tang F, Huang Y. Single-cell RNA-seq transcriptome analysis of linear and circular RNAs in mouse preimplantation embryos. Genome Biol. 2015;16(1):148.

    PubMed  PubMed Central  Google Scholar 

  30. Szabo L, Morey R, Palpant NJ, Wang PL, Afari N, Jiang C, Parast MM, Murry CE, Laurent LC, Salzman J. Statistically based splicing detection reveals neural enrichment and tissue-specific induction of circular RNA during human fetal development. Genome Biol. 2015;16(1):126.

    PubMed  PubMed Central  Google Scholar 

  31. Verma N, Rettenmeier AW, Schmitz-Spanke S. Recent advances in the use of Sus scrofa (pig) as a model system for proteomic studies. Proteomics. 2011;11(4):776–93.

    CAS  PubMed  Google Scholar 

  32. Cen J, Liang Y, Huang Y, Pan Y, Shu G, Zheng Z, Liao X, Zhou M, Chen D, Fang Y, et al. Circular RNA circSDHC serves as a sponge for miR–127–3p to promote the proliferation and metastasis of renal cell carcinoma via the CDKN3/E2F1 axis. Mol Cancer. 2021;20(1):19.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.

    CAS  PubMed  Google Scholar 

  34. Qi K, Liu Y, Li C, Li X, Li X, Wang K, Qiao R, Han X. Construction of circRNA-related ceRNA networks in longissimus dorsi muscle of Queshan Black and large White pigs. Mol Genet Genomics. 2022;297(1):101–12.

    CAS  PubMed  Google Scholar 

  35. Li M, Zhang N, Li J, Ji M, Zhao T, An J, Cai C, Yang Y, Gao P, Cao G, et al. CircRNA profiling of skeletal muscle in two Pig Breeds reveals CircIGF1R regulates myoblast differentiation via miR–16. Int J Mol Sci. 2023;24(4):3779.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. Cesta MF. Normal structure, function, and histology of the spleen. Toxicol Pathol. 2006;34(5):455–65.

    PubMed  Google Scholar 

  37. Dardalhon V, Schubart AS, Reddy J, Meyers JH, Monney L, Sabatos CA, Ahuja R, Nguyen K, Freeman GJ, Greenfield EA, et al. CD226 is specifically expressed on the surface of Th1 cells and regulates their expansion and effector functions. J Immunol. 2005;175(3):1558–65.

    CAS  PubMed  Google Scholar 

  38. Horejsí V, Zhang W, Schraven B. Transmembrane adaptor proteins: organizers of immunoreceptor signalling. Nat Rev Immunol. 2004;4(8):603–16.

    PubMed  Google Scholar 

  39. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA–MEM. arXiv. 2013;26.

  40. Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):4.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    CAS  PubMed  Google Scholar 

  42. Jin Y, Hammell M. Analysis of RNA-Seq Data using TEtranscripts. Methods Mol Biol. 2018;1751:153–67.

    CAS  PubMed  Google Scholar 

  43. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    PubMed  PubMed Central  Google Scholar 

  44. Smedley D, Haider S, Ballester B, et al. BioMart–biological queries made easy. BMC Genomics. 2009;10:22.

    PubMed  PubMed Central  Google Scholar 

  45. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.

    PubMed  PubMed Central  Google Scholar 

  46. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 202;51:D587–92.

  48. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2(3):100141.

    CAS  Google Scholar 

Download references

Acknowledgements

We thank Oebiotech Corporation (Shanghai, China) for Illumina sequencing.

Funding

This study was supported by the Seed Industry Vitalization Research Projects of Jiangsu Province (JBGS[2021]099); Jiangsu Agricultural Science and Technology Innovation Fund (SCX(23)3143); the Open Project of International Research Laboratory of Prevention and Control of Important Animal Infectious Diseases and Zoonotic Diseases of Jiangsu Higher Education Institutions, Yangzhou University, Yangzhou, China; China Postdoctoral Science Foundation (2021M702764) and the Priority Academic Program Development of Jiangsu Higher Education Institutions.

Author information

Author notes

  1. Yifu Wang and Jinhua Cheng contributed equally to this work.

Authors and Affiliations

  1. College of Animal Science and Technology, Yangzhou University, Yangzhou, 225009, China

    Yifu Wang, Chao Xu, Wenbin Bao, Shenglong Wu & Zhengchang Wu

  2. Institute of Comparative Medicine, College of Veterinary Medicine, Yangzhou University, Yangzhou, 225009, China

    Jian Jin & Zhengchang Wu

  3. International Research Laboratory of Prevention and Control of Important Animal Infectious Diseases and Zoonotic Diseases of Jiangsu Higher Education Institutions, Yangzhou University, Yangzhou, 225009, China

    Zhengchang Wu

  4. Institute of Animal Science, Jiangsu Academy of Agricultural Sciences, Nanjing, 210014, China

    Jinhua Cheng

Authors
  1. Yifu Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jinhua Cheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Chao Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jian Jin
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Wenbin Bao
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Shenglong Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Zhengchang Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Yifu Wang: Conceptualization, Writing-original draft, Validation. Jinhua Cheng: Resources, Methodology, Validation. Chao Xu: Data curation, Software. Jian Jin: Validation. Wenbin Bao: Project administration, Writing-review & editing. Shenglong Wu: Project administration. Zhengchang Wu: Supervision, Validation, Writing-original draft. All authors reviewed the manuscript.

Corresponding author

Correspondence to Zhengchang Wu.

Ethics declarations

Ethics approval and consent to participate

All experiments were approved by the Institutional Animal Care and Use Committee (IACUC) of Yangzhou University (Pig: SYXK(Su)2021-0026) and were performed according to the Animal Ethics Procedures and Guidelines of the People’s Republic of China. We confirm that the study was conducted in accordance with ARRIVE guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, Y., Cheng, J., Xu, C. et al. Genome-wide analysis of circRNA regulation during spleen development of Chinese indigenous breed Meishan pigs. BMC Genomics 24, 477 (2023). https://doi.org/10.1186/s12864-023-09612-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09612-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us