Comprehensive Analysis of CircRNA Expression Profiles in Multiple Tissues of Pigs

Circular RNAs (circRNAs) are a class of non-coding RNAs with diverse functions, and previous studies have reported that circRNAs are involved in the growth and development of pigs. However, studies about porcine circRNAs over the past few years have focused on a limited number of tissues. Based on 215 publicly available RNA sequencing (RNA-seq) samples, we conducted a comprehensive analysis of circRNAs in nine pig tissues, namely, the gallbladder, heart, liver, longissimus dorsi, lung, ovary, pituitary, skeletal muscle, and spleen. Here, we identified a total of 82,528 circRNAs and discovered 3818 novel circRNAs that were not reported in the CircAtlas database. Moreover, we obtained 492 housekeeping circRNAs and 3489 tissue-specific circRNAs. The housekeeping circRNAs were enriched in signaling pathways regulating basic biological tissue activities, such as chromatin remodeling, nuclear-transcribed mRNA catabolic process, and protein methylation. The tissue-specific circRNAs were enriched in signaling pathways related to tissue-specific functions, such as muscle system process in skeletal muscle, cilium organization in pituitary, and cortical cytoskeleton in ovary. Through weighted gene co-expression network analysis, we identified 14 modules comprising 1377 hub circRNAs. Additionally, we explored circRNA–miRNA–mRNA networks to elucidate the interaction relationships between tissue-specific circRNAs and tissue-specific genes. Furthermore, our conservation analysis revealed that 19.29% of circRNAs in pigs shared homologous positions with their counterparts in humans. In summary, this extensive profiling of housekeeping, tissue-specific, and co-expressed circRNAs provides valuable insights into understanding the molecular mechanisms of pig transcriptional expression, ultimately deepening our understanding of genetic and biological processes.


Introduction
Circular RNA (circRNA) is a type of non-coding RNA generated through a unique splicing process known as backsplicing [1].Backsplicing is a specific mechanism that produces downstream splice-donor sites and upstream splice-acceptor sites, resulting in their covalent connection [2].Various types of circRNAs have been identified, including exon, intron, and exon-intron type [3].Most circRNAs (exonic circRNAs) are expressed from known protein-coding genes and consist of a single exon or of multiple exons [4,5].Additionally, a minority of circRNAs (intronic circRNAs) are composed of introns [6].Other circRNAs (exon-intron circRNAs) contain sequences derived from both exons and introns due to internal intron retention or a failure in the debranching of intronic lariats during canonical splicing [7][8][9].The circular structure of circRNA gives it higher stability and resistance to degradation.This particular structure of circRNA allows circRNAs to persist and remain stable in cells for extended periods of time [10].Meanwhile, circRNAs have a range of biological functions, including competitive miRNA inhibition, direct regulation of parental mRNA expression, and protein sequestration [11][12][13].
Pigs exhibit a significant degree of physiology and anatomy similarity to humans, rendering them a favorable animal model for investigating human diseases and developmental processes [14].Moreover, domestic pigs are essential farm animals with a long history of selective breeding and immense economic value [15].Over time, humans have conducted extensive breeding programs for domestic pigs, resulting in a rich genetic diversity and the development of various desirable traits [16].Research on circRNAs is an important area of study for pig breeding, as differential expression of circRNAs has been shown to impact various traits, including meat quality, reproductive ability, and disease susceptibility.For instance, Li et al. found that circIGF1R could act as a miRNA sponge to negatively regulate miR-16 to promote the myoblast differentiation of the porcine skeletal muscle satellite cells, thereby influencing meat quality [17].Wang et al. detected 60,78, and 86 differentially expressed circRNAs that played a key role in muscle development and lipid deposition across 38, 58, and 78 days post conception [18].Zhuang et al. identified circKANSL1L, which influences muscle fiber type differentiation and exhibits a high level of conservation between mice and pigs [19].Niu et al. detected 305 differentially expressed circRNA in porcine ovaries, and they were enriched in several reproductive-related signaling pathways [20].Mester-Tonczar et al. demonstrated that circRNA-CDR1as beneficially impacts cardiac function in pigs by down-regulating miR-7 in the heart [21].Although there have been some advances in pig circRNA research over the past few years, investigations are still mainly focused on a limited number of tissues.
The effect of circRNA on pig traits constitutes a complex regulatory network involving multiple processes and mechanisms.A better comprehension of the circRNA profiles in different pig tissues will help us to leverage the pig genome for applications in various fields, such as species of evolution, growth and development, breeding, construction of models for humans, and gene-disease phenotypic association analysis.With the advancement of RNA sequencing technology, an increasing number of pig genomic data have been made publicly available and widely utilized [22].Numerous studies of non-coding RNA based on RNA-seq and the related databases have been published [23,24].The databases of circRNA were published as well, such as CircFunBase and CircAtlas [25,26], but the number of circRNAs in pigs should be greatly enriched.Furthermore, recent studies have found a significant amount of circRNA in porcine muti-tissues.For example, Jin et al. found 48,232 circRNAs from 31 tissues and two cell lines [27].However, circRNA expression in multiple tissues of pigs should still be explored based on bulk data.
In this study, we aimed to investigate the circRNA profiles by RNA-seq from nine porcine tissues (gallbladder, heart, liver, longissimus dorsi, lung, ovary, pituitary, skeletal muscle, and spleen) to explore the association and distinctness of circRNA expression patterns among different tissues.Our observations will provide new insights into the biological processes of various tissues in pigs.

CircRNA Detection in Multiple Tissues
CircRNAs were detected using the CIRI2 and find_cir methods in a total of 215 samples from nine different tissues (see Section 2 for details).In summary, 238,923 circRNAs were detected using the find_cir, while 111,619 circRNAs were identified utilizing the CIRI2.Among these, a total of 82,528 circRNAs were detected by both methods and used for subsequent analysis (Figure 1a, Table S2).Compared with the CircAtlas database (https://ngdc.cncb.ac.cn/circatlas/, accessed on 7 November 2023) [26], in total, 3818 new circRNA were found (Table S3).The tissue-specific distribution of circRNAs indicated that the highest number (14,224) of circRNAs was detected in the pituitary, while the lowest number (1709) was detected in the heart compared to other tissues (Figure 1b).Additionally, the longissimus dorsi exhibited the highest mean junction ratio (0.15) among the tissues, while the spleen showed the lowest mean junction ratio (0.06) compared to other tissues (Figure 1c).The average genome distance of all circRNAs was calculated to be 15,383 bp, with circRNAs shorter than 50,000 bp accounting for 93.9% of the total (Figure 1d).Notably, the majority of circRNAs belonged to the exonic type, accounting for 88.0% (pituitary) to 92.6% (gallbladder) across different tissues.The remaining circRNA types, including intronic and intergenic, accounted for 4.2% (gallbladder) to 7.4% (skeletal muscle) and 3.2% (gallbladder) to 6.8% (heart), respectively (Figure 1e).
circRNAs was detected in the pituitary, while the lowest number (1709) was detected in the heart compared to other tissues (Figure 1b).Additionally, the longissimus dorsi exhibited the highest mean junction ratio (0.15) among the tissues, while the spleen showed the lowest mean junction ratio (0.06) compared to other tissues (Figure 1c).The average genome distance of all circRNAs was calculated to be 15,383 bp, with circRNAs shorter than 50,000 bp accounting for 93.9% of the total (Figure 1d).Notably, the majority of circRNAs belonged to the exonic type, accounting for 88.0% (pituitary) to 92.6% (gallbladder) across different tissues.The remaining circRNA types, including intronic and intergenic, accounted for 4.2% (gallbladder) to 7.4% (skeletal muscle) and 3.2% (gallbladder) to 6.8% (heart), respectively (Figure 1e).

CircRNA Clustering Analysis
Clustering analysis was performed on the retained circRNAs, and the clustering results demonstrated that most tissues exhibited well-defined clusters.However, a noticeable stratification was observed within the skeletal muscle tissues (Figure 2a).Principal component analysis (PCA) showed that the skeletal muscle samples could be classified into two distinct groups based on their developmental stages: embryonic muscle and postnatal muscle (Figure S1).
Furthermore, circRNA average expression levels across tissues were used to perform Spearman clustering analysis.Specifically, a strong correlation was observed between skeletal muscle and the longissimus dorsi (r = 0.55), indicating a similarity in the expression profiles of circRNAs in these tissues.Similarly, a significant correlation was found between the ovary and pituitary gland (r = 0.53) (Figure 2b).
Through an analysis of the parental genes associated with circRNAs between tissues, it was observed that more than 99% (1427/1438) of circRNAs in the longissimus dorsi tissue were shared with skeletal muscle (Figure 2c).On the other hand, circRNAs shared by

CircRNA Clustering Analysis
Clustering analysis was performed on the retained circRNAs, and the clustering results demonstrated that most tissues exhibited well-defined clusters.However, a noticeable stratification was observed within the skeletal muscle tissues (Figure 2a).Principal component analysis (PCA) showed that the skeletal muscle samples could be classified into two distinct groups based on their developmental stages: embryonic muscle and postnatal muscle (Figure S1).
Furthermore, circRNA average expression levels across tissues were used to perform Spearman clustering analysis.Specifically, a strong correlation was observed between skeletal muscle and the longissimus dorsi (r = 0.55), indicating a similarity in the expression profiles of circRNAs in these tissues.Similarly, a significant correlation was found between the ovary and pituitary gland (r = 0.53) (Figure 2b).
Through an analysis of the parental genes associated with circRNAs between tissues, it was observed that more than 99% (1427/1438) of circRNAs in the longissimus dorsi tissue were shared with skeletal muscle (Figure 2c).On the other hand, circRNAs shared by the pituitary gland and ovary exhibited a high degree of overlap in their parental genes, with over 50% (2586/4091) of these host genes being shared with both tissues (Figure 2d).
the pituitary gland and ovary exhibited a high degree of overlap in their parental genes, with over 50% (2586/4091) of these host genes being shared with both tissues (Figure 2d).

Tissue-Specific and Housekeeping CircRNA Profile
After a strict screening (see Section 2 for detail), in total, 24,745 circRNAs were identified by pervious materials in this study, of which 492 were categorized as housekeeping circRNAs and 3489 were classified as tissue-specific circRNAs (Figure 3a, Table S4).Furthermore, 403 of the 492 housekeeping circRNAs' parental genes were housekeeping genes identified by the pigGtex project (Table S5).To gain insights into the biological functions of these circRNAs, the parental genes (housekeeping circRNAs) associated with them were subjected to GO pathway analyses.The GO enrichment analysis of housekeeping circRNAs enriched in chromatin remodeling and nuclear-transcribed mRNA catabolic processes in the biological process category, like histone binding (Figure 3b, Table S6), suggested that housekeeping circRNAs may be involved in gene expression regulation and RNA metabolism.In the tissue-specific circRNAs, the parental genes are also able to be enriched for their associated functions, for instance, muscle system process in skeletal muscle and neuron-to-neuron synapse in pituitary (Figure S2, Table S7).
We found that the tissue-specific circRNAs were significantly related to traits of meat and carcass, and health and reproduction (Table S10).Additionally, the KEGG pathway analysis of housekeeping circRNAs revealed significant enrichment, such as long-term depression, long-term potentiation, TGF-beta signaling pathway, and autophagy (Figure 3c, Table S8).These pathways are known to play crucial roles in various cellular processes and signaling cascades.Similarly, the parent genes of tissue-specific circRNA are also able to be enriched for their associated functions (Figure S3).
Furthermore, the miRNAs with the top five highest scores in miRanda-based circR-NAs match were selected as potential miRNA targets (see Section 2 for detail).We found

Tissue-Specific and Housekeeping CircRNA Profile
After a strict screening (see Section 2 for detail), in total, 24,745 circRNAs were identified by pervious materials in this study, of which 492 were categorized as housekeeping circRNAs and 3489 were classified as tissue-specific circRNAs (Figure 3a, Table S4).Furthermore, 403 of the 492 housekeeping circRNAs' parental genes were housekeeping genes identified by the pigGtex project (Table S5).To gain insights into the biological functions of these circRNAs, the parental genes (housekeeping circRNAs) associated with them were subjected to GO pathway analyses.The GO enrichment analysis of housekeeping circRNAs enriched in chromatin remodeling and nuclear-transcribed mRNA catabolic processes in the biological process category, like histone binding (Figure 3b, Table S6), suggested that housekeeping circRNAs may be involved in gene expression regulation and RNA metabolism.In the tissue-specific circRNAs, the parental genes are also able to be enriched for their associated functions, for instance, muscle system process in skeletal muscle and neuron-to-neuron synapse in pituitary (Figure S2, Table S7).
We found that the tissue-specific circRNAs were significantly related to traits of meat and carcass, and health and reproduction (Table S10).Additionally, the KEGG pathway analysis of housekeeping circRNAs revealed significant enrichment, such as long-term depression, long-term potentiation, TGF-beta signaling pathway, and autophagy (Figure 3c, Table S8).These pathways are known to play crucial roles in various cellular processes and signaling cascades.Similarly, the parent genes of tissue-specific circRNA are also able to be enriched for their associated functions (Figure S3).
Furthermore, the miRNAs with the top five highest scores in miRanda-based circRNAs match were selected as potential miRNA targets (see Section 2 for detail).We found that tissue-specific circRNA might interact with many miRNAs or indirectly with differentially expressed tissue-specific genes in liver, longissimus dorsi, lung, ovary, and pituitary tissues (Figure S4).For instance, we found that "circ X:27122083-27142550" might interact with that tissue-specific circRNA might interact with many miRNAs or indirectly with differentially expressed tissue-specific genes in liver, longissimus dorsi, lung, ovary, and pituitary tissues (Figure S4).For instance, we found that "circ X:27122083-27142550" might interact with STAT2 via ssc-mir-615, might interact with CCL4 via ssc-mir-1842, and might interact with FETUB via ssc-mir-7-2 in liver (Figure 3d).

Co-Expression CircRNA Network across Pig Tissues
WGCNA was performed on the filtered circRNAs to investigate the biological relationships and potential functions of the core circRNAs among tissues.In this analysis, a soft threshold value (β) was determined to construct an efficient scale-free network.When the R 2 (coefficient of determination) value exceeded 0.8, the threshold value of 4 was selected (Figure S5).By applying the hierarchical clustering algorithm, in total, 16 distinct modules were obtained.These modules were subsequently corrected and merged, resulting in 14 modules (Figures 4a and S5).Each module represents a set of circRNAs that exhibit highly correlated expression patterns across tissues.Among these modules, seven module-tissue relationships were tissue-specific modules (r > 0.65).For instance, the salmon module displayed a strong correlation with heart (r = 0.96), and the red module showed a strong correlation with liver (r = 0.93), and the magenta exhibited a strong correlation with lung (r = 0.89) (Figure 4b).Furthermore, among these seven tissue-specific modules, we applied the MM and GS methods to identify hub-circRNA genes by using a cutoff of |GS| > 0.2 and |MM| > 0.8.Finally, 142, 262, 160, 149, 29, 625, and 10 hub circRNAs were observed in heart-salmon, liver-red, lung-magenta, skeletal muscle-black, skeletal muscle-green, pituitaryturquoise, and spleen-blue modules, respectively (Table S9).

Co-Expression CircRNA Network across Pig Tissues
WGCNA was performed on the filtered circRNAs to investigate the biological relationships and potential functions of the core circRNAs among tissues.In this analysis, a soft threshold value (β) was determined to construct an efficient scale-free network.When the R 2 (coefficient of determination) value exceeded 0.8, the threshold value of 4 was selected (Figure S5).By applying the hierarchical clustering algorithm, in total, 16 distinct modules were obtained.These modules were subsequently corrected and merged, resulting in 14 modules (Figure 4a and Figure S5).Each module represents a set of circRNAs that exhibit highly correlated expression patterns across tissues.Among these modules, seven module-tissue relationships were tissue-specific modules (r > 0.65).For instance, the salmon module displayed a strong correlation with heart (r = 0.96), and the red module showed a strong correlation with liver (r = 0.93), and the magenta exhibited a strong correlation with lung (r = 0.89) (Figure 4b).

Conservation between Pig and Human
We compared the circRNAs between the pig and human genomes using the LiftOver tool (minMatch = 0.5).We found that 77.10% of circRNAs (19,078/24,745) in the pig genome could be aligned with sequences in the human genome.Furthermore, 19.29% of the pig circRNAs (4774/24,745) were sequenced and usage conserved with circRNAs in the human genome (Table 1).These conserved circRNAs are all derived from orthologous genes between pig and human (Table S11), and 96% of the conserved circRNA originate from exons (Figure S6).

Discussion
Pigs serve as an essential economic resource and an ideal animal model for studying various aspects of animal domestication and human diseases.In recent years, circRNAs have emerged as a prevalent class of non-coding RNAs with diverse functions.However, the pattern of circRNA expression based on large-scale data from diverse tissues remains unclear in pigs.In our study, we aimed to comprehensively explore the landscape of circR-NAs in pigs by analyzing a dataset comprising nine different pig tissues.In total, 3489 tissue-specific circRNAs and 492 housekeeping circRNAs were preliminarily identified.Meanwhile, 14 modules with 1377 hub circRNAs were obtained.
To ensure the accuracy of the results, we employed two well-established algorithms, CIRI2 and Find_cir, to detect circRNAs using RNA sequencing data.CIRI2 algorithm leverages paired chiastic clipping signals obtained from the mapping information, while Find_cir predicts back-splicing events by examining the first and last 20 bp anchors of unmapped reads [31,32].These algorithms have been extensively validated and used in previous studies [50].It is important to note that due to differences in their underlying principles and strategies, these algorithms may yield different predictions for circRNAs.

Conservation between Pig and Human
We compared the circRNAs between the pig and human genomes using the LiftOver tool (minMatch = 0.5).We found that 77.10% of circRNAs (19,078/24,745) in the pig genome could be aligned with sequences in the human genome.Furthermore, 19.29% of the pig circRNAs (4774/24,745) were sequenced and usage conserved with circRNAs in the human genome (Table 1).These conserved circRNAs are all derived from orthologous genes between pig and human (Table S11), and 96% of the conserved circRNA originate from exons (Figure S6).

Discussion
Pigs serve as an essential economic resource and an ideal animal model for studying various aspects of animal domestication and human diseases.In recent years, circRNAs have emerged as a prevalent class of non-coding RNAs with diverse functions.However, the pattern of circRNA expression based on large-scale data from diverse tissues remains unclear in pigs.In our study, we aimed to comprehensively explore the landscape of circRNAs in pigs by analyzing a dataset comprising nine different pig tissues.In total, 3489 tissue-specific circRNAs and 492 housekeeping circRNAs were preliminarily identified.Meanwhile, 14 modules with 1377 hub circRNAs were obtained.
To ensure the accuracy of the results, we employed two well-established algorithms, CIRI2 and Find_cir, to detect circRNAs using RNA sequencing data.CIRI2 algorithm leverages paired chiastic clipping signals obtained from the mapping information, while Find_cir predicts back-splicing events by examining the first and last 20 bp anchors of unmapped reads [28,29].These algorithms have been extensively validated and used in previous studies [30].It is important to note that due to differences in their underlying principles and strategies, these algorithms may yield different predictions for circRNAs.Hence, we adopted a conservative approach by combining the outputs of both algorithms to minimize false-positive predictions and enhance the reliability of our circRNA dataset.By this approach, more than 95% of detected circRNA could be found in the released database which could certify the accuracy of this method [31].
Statistical analysis of the detected circRNAs revealed that exon-derived circRNAs accounted for approximately 90% of the total circRNAs across pig tissues.This observation aligns with previous studies in pigs and other species [30,32].Additionally, we observed that the distribution of circRNAs peaked within a genomic span of 2000 bp, with a gradual decrease in abundance as the genomic span increased.This distribution pattern is consistent with known characteristics of circRNAs [30,33].Furthermore, our study revealed tissuespecific differences in the quantity and splicing rates of circRNAs.This suggests that circRNA biogenesis and regulation are subject to tissue-specific factors, such as alternative splicing and post-transcriptional processing.These tissue-specific differences highlight the potential functional diversity and regulatory roles of circRNAs in different biological contexts [34,35].
The expression analysis of circRNAs in pigs revealed distinct tissue-specific expression profiles across different tissues.Subsequently, the study further investigated and identified 3489 tissue-specific circRNAs and 492 housekeeping circRNAs.By performing pathway enrichment analysis on the parental genes of tissue-specific circRNAs, key pathways related to tissue-specific functions were uncovered.Taking the heart tissue as an example, several parental genes were found to be involved in the formation of intercellular junction complexes.For instance, the parental genes PKP2 of "circ5:41492479|41493703" and "circ5:41438437|41447707" and CDH2 of "circ6:112408292|112462754" interact with other proteins, such as intercellular adhesion proteins and cytoskeletal proteins, enhancing the adhesion and stability between cardiac muscle cells [36,37].In the liver, the parental gene GYS2 of "circ5:51884381|51889554" plays a role in glycogen synthesis and metabolic regulation [38].In the longissimus dorsi, the parental gene CACNA1S of "circ10:23546471|23549237" encodes the α1S subunit of a voltage-gated calcium ion channel, which is predominantly expressed in skeletal muscle cells and regulates muscle contraction [39].In the lungs, the parental gene CTSS of "circ4:98422323|98425826" is associated with airway inflammation and lung disease development, including its involvement in lung cancer progression [40].In the ovaries, the parental gene PBX1 of "circ4:85892566|85906932" plays a critical role in embryonic development, contributing to processes such as axon guidance, organ formation, and segmental differentiation [41].In the pituitary gland, the parental gene ADCY9 of "circ3:38275930|38288061" encodes adenylate cyclase 9, which catalyzes the cyclization of adenosine monophosphate and generates the second messenger cyclic adenosine monophosphate (cAMP) [42].In skeletal muscle, the parental gene ACTN2 of "circ14:54704317|54716478" encodes alpha-actinin-2, a myofibrillar protein responsible for constructing and maintaining the cell cytoskeleton.ACTN2 plays a crucial role in muscle contraction and force transmission [43].In the spleen, the parental gene ITGAL of "circ3:17840569|17841414" regulates the activity and function of immune cells by participating in their adhesion, migration, and interactions within the spleen [44].Notably, specific circRNAs in the gallbladder tissue, such as SLC4A4 and SLC5A1, are primarily involved in gastrointestinal tract-related functions [45,46].These findings highlight the functional diversity and tissue-specific roles of circRNAs in different organs.The identified parental genes associated with tissue-specific circRNAs provide insights into the potential regulatory mechanisms and biological functions of circRNAs in specific tissues.
Similarly, the analysis of housekeeping circRNAs and their parental genes revealed several circular RNAs that are involved in key pathways related to fundamental cellular functions.For instance, the parental gene MAP2K1 of "circ1:164420445|164422619" encodes mitogen-activated protein kinase 1 (MAP2K1), which participates in the activation of the MAPK signaling pathway [47].This pathway regulates critical biological processes such as cell growth, differentiation, proliferation, and survival.Another example is the parental gene ATG4C of "circ6:149633393|149666969", which is involved in the regulation of autophagy, a highly regulated process essential for maintaining cellular survival and metabolic balance [48].Autophagy plays a crucial role in cellular homeostasis by recycling damaged organelles and proteins.Additionally, the parental gene CREBBP of "circ3:38413987|38414699" and "circ3:38413987|38450234" is involved in the regulation of the cell cycle progression.CREBBP controls cell proliferation and differentiation fate, and it can interact with various apoptosis-regulating proteins, modulating the apoptotic signaling pathway in cells [49].These findings indicate that housekeeping circRNAs derived from these parental genes are involved in important cellular pathways.These circRNAs likely contribute to the regulation of essential cellular functions, including cell growth, survival, metabolism, and differentiation.The identification of these housekeeping circRNAs and their associated parental genes sheds light on the functional roles of circRNAs in fundamental cellular processes.Understanding the regulatory mechanisms and functions of these housekeeping circRNAs can provide valuable insights into the molecular mechanisms underlying cellular homeostasis and cellular processes associated with development, differentiation, and disease.
Additionally, we investigated the binding capabilities of circRNAs to miRNAs based on their sequence complementarity.By evaluating the binding affinity between miRNAs and circRNAs, we identified a subset of circRNAs with strong binding capabilities, suggesting their potential role as competitive endogenous RNAs (ceRNAs).CeRNAs can sequester miRNAs and influence mRNA expression, forming a complex regulatory ceRNA network [19].For example, in the liver, ebv-circLMP2A exists by regulating miR-3908, and miR-3908 further modulates the specifically expressed STAT2 [50,51].
Furthermore, we performed co-expression circRNA network analysis to examine the interactions among circRNAs.The analysis revealed that all circRNAs could be classified into 14 distinct modules, each representing a set of co-expressed circRNAs.Notably, nine of these modules exhibited strong correlations with individual tissues, suggesting their involvement in tissue-specific functions.From these modules, we identified hub circRNAs that likely play crucial roles in tissue-specific functions.Hub circRNAs are highly connected within the co-expression network and are considered key regulators within the ceRNA network.The identification of hub circRNAs provides valuable insights into the potential regulatory roles of circRNAs in tissue-specific processes.
Moreover, we conducted a conservation analysis to assess the conservation of cir-cRNAs between humans and pigs.The results indicated that approximately 20% of the identified circRNAs share homologous genomic positions between the two species.The percentage is consistent with conservation studies between human and mouse [52].This conservation suggests that these circRNAs may possess similar sequence and functional characteristics across species, highlighting their potential importance in cross-species biological processes.
Finally, although 24,745 circRNAs were identified in this study, the data volume needs further supplementation.There are a large number of genomic data in pigGTEx, but circRNA research is still lacking [53].In the future, we hope to include more data, such as pigGTEx data, in the circRNA analysis, aiming to discover more novel circRNAs and thus more comprehensively complement the circRNA data.
Overall, our findings provide valuable insights into the circRNA expression profiles and their potential roles in the biological development and processes of various pig tissues.This information lays the foundation for further investigations into the functional significance and regulatory mechanisms of circRNAs in pigs, ultimately contributing to our understanding of tissue-specific gene regulation and the molecular basis of pig biology.
However, there are still some limitations in this study.Although Teng et al. have demonstrated that tissue-specific factors have a greater influence in transcriptome analysis than sex, breeding, age, and other factors [53], we categorized these factors as batch effects in our study due to the absence of some annotation information in the data.This approach, as compared to individually correcting for each factor, may introduce some deviation.Furthermore, the lack of corresponding phenotypic data means that we can only indirectly assess the relationship between circRNAs and phenotypes through association analysis between circRNAs and known trait-related QTLs.
In future research, we hope to have more comprehensive annotated data to obtain a more comprehensive pig circRNA profile, allowing for a more in-depth investigation into the relationship between circRNA expression and phenotypes.S1.The pig reference genome (Sscrofa11.1) was used for the analysis.

Quality Control and Alignment
The samples were transformed into pair-end fastq format using the fast-dump module in the sra-toolkit (version 2.8.2), and the quality of the raw sequencing data was assessed using FASTP software (version 0.23.2) with default parameters to filter and trim the reads to remove low-quality bases and adaptors [54].Subsequently, the clean data were mapped to the reference genome using HISAT2 software (version 2.1.0),and the SAM files were converted to BAM format using Samtools (version v1.11) [55,56].

CircRNA Detection
To ensure the most reliable precision and sensitivity in circRNA identification, a combination of CIRI2 and find_cir methods were employed on the clean data [28,29].In order to eliminate the impact of variables such as sex, breeding, age, and sequencing methods, we integrated them into batch effects and subsequently applied the combat algorithm from the R package sva after grouping the samples by their source for batch effect correction [57].Subsequently, the number of circRNAs, their length, and the junction ratio (the proportion of circular junction reads compared to linear junction reads) were counted and quantified for each tissue.Additionally, three circRNA types were detected in each tissue, including exonic circRNAs originating from exons, intronic circRNAs originating from introns, and intergenetic circRNA originating from the intergenetic region.

CircRNA Profile between Tissues
To minimize potentially spurious events, only circRNAs detected in more than 50% of the samples within at least one tissue were considered for further analysis [58].This filtering criterion resulted in a total of 24,745 circRNAs meeting the inclusion criteria.Subsequently, a t-distributed stochastic neighbor embedding (t-SNE) analysis was performed to explore the relationship and patterns among the samples [59], and based on the expression levels of circRNAs to assess the correlation between tissues, the Spearman's correlation coefficient was calculated [60].

Identification of Housekeeping and Tissue-Specific CircRNAs
Two categories were defined to categorize circRNAs based on their expression patterns, housekeeping circRNAs, and tissue-specific circRNAs.The housekeeping circRNAs were identified by the conditions for the housekeeping mRNA identification, as in a previous study [61].On the other hand, tissue-specific circRNAs were defined as those expressed exclusively in one tissue.To visualize the expression patterns of circRNAs across tissues, a heatmap was generated using the Pheatmap package of R (version 1.0.12).

GO and KEGG Enrichment
Two tools were utilized to analyze the parental genes of the identified circRNAs and gain insights into their potential functional roles, the clusterProfiler package of R for Gene Ontology (GO) enrichment analysis and the ShinyGO online server (http:// bioinformatics.sdstate.edu/go/accessed on 1 January 2023) for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis [62,63].In both analyses, the hypergeometric test was applied to determine the significance of enrichment.A significance threshold of p < 0.05 was considered, indicating that genes associated with circRNAs were significantly enriched in specific GO terms or KEGG pathways.

Competing Endogenous RNA Network Construction
To investigate the potential regulatory interactions between circRNAs (housekeeping circRNAs and tissue-specific circRNAs), miRNAs, and mRNAs, miRanda software (version 3.3a) was used to predict the interaction between circRNA-miRNA and miRNA-mRNA, a match score above 140 was used as a threshold to predict significant interactions [64].Firstly, the miRNA sequences were obtained from the miRBase database (http://www.mirbase.org/accessed on 1 March 2023) [65], and the sequences of cir-cRNAs and mRNAs were obtained from a previous study using bedtools2 software (version 2.23.0) [66,67].Subsequently, miRanda was used to calculate the free energy and score between miRNAs and circRNAs, as well as between miRNAs and mRNAs.These calculations provided information about the stability and potential strength of interactions between miRNAs and their targets.Finally, the regulatory interactions among circRNA-miRNA-mRNA were visualized using Cytoscape software (version 3.10.0)[68].

Co-Expression Network Analysis
The Weighted Gene Co-expression Network Analysis (WGCNA) method was employed to investigate the co-expression patterns of circRNAs among tissues [69].The hclust function with the average agglomeration method for the WGCNA package was used for cluster analysis.The pickSoftThreshold function from the WGCNA package was utilized to determine the soft threshold (β), which determines the connectivity strength between circRNAs.The module with a correlation coefficient larger than 0.65 was considered the important tissue-specific module in this study [61].Subsequently, the important tissuespecific modules were identified from the co-expression network and further analyzed to identify core circRNAs within them.Briefly, the key circRNAs within the tissue-specific modules, which are called hub circRNAs, were identified using Module Membership (MM) and Gene Significance (GS) methods.Hub circRNAs were determined based on predefined thresholds, typically |GS| > 0.2 and |MM| > 0.8, indicating strong correlations between circRNAs and module eigengenes [70].

Association Analysis of Traits
To identify the relationship between tissue-specific circRNAs and phenotypes, the QTL data of pigs were downloaded from the QTLdb database (https://www.animalgenome.org/cgi-bin/QTLdb/index accessed on 1 November 2023), and we counted any QTLs that overlapped with circRNAs in terms of location and subsequently performed Fisher's tests to obtain the degree of association between circRNA and trait-related QTL [71].

Conservation Analysis
To identify and compare circRNAs between pig and human, the predicted circR-NAs were converted from pig genomic locations (susScr11) to human genomic locations (hg19) by using the LiftOver tool from the UCSC genome browser with the parameter "minMatch = 0.5" [72].The circRNAs aligned to the human genome were then compared with the known circRNAs in humans, which were obtained from the circBase database (http://www.circbase.org/accessed on 1 May 2023) [31].

Conclusions
In summary, our study focused on exploring the expression profiles of circRNAs in nine tissues of pigs.We successfully identified 492 housekeeping circRNAs and 3489 tissue-specific circRNAs, highlighting their distinct expression patterns and potential roles in tissue-specific processes.Moreover, we obtained 14 modules consisting of 1210 hub circRNAs which are likely to play crucial regulatory roles in various biological activities across tissues.Notably, we have discovered a total of 3818 novel circRNAs, contributing additional information to augment the pig circRNA dataset.

Table 1 .
Number of conservation circRNAs between pig and human.

Table 1 .
Number of conservation circRNAs between pig and human.