Content area
Background
Circular RNAs (circRNAs) and mRNAs are distinct transcripts from the same genes, produced by different splicing mechanisms. This study investigates the behavior of the circular transcriptome relative to the linear one across biological conditions and tissues. We analyzed transcriptomic data from 36 bovine monocyte-derived macrophage (MDM) samples collected during an ex vivo Mycobacterium avium ssp. paratuberculosis (MAP) infection experiment, stratified by Johne’s disease (JD) antibody status (JD+ or JD−) and by infection condition (control or MAP infected). We extended our analysis to healthy bovine tissues, including neonatal and post-pubertal testes, and liver and muscle samples from 12 animals stratified by sex and feed efficiency.
Results
In the 36 MDM samples, we identified 3358 exonic circRNAs derived from 1895 genes. By comparing the mean expression levels of circRNAs and linear transcripts, and considering the number of expressed genes, we estimate that the circular transcriptome is approximately 100 times smaller than the linear transcriptome. Analyses of the circular and linear transcriptomes revealed that MAP infection impacted only the linear transcriptome of MDM_JD− . The other three transcriptomes—circular JD− , circular JD+ , and linear JD+ —showed no infection-specific response. In the testes, maturation was associated with profound but uncoordinated changes in the circular and linear transcriptomes. While circRNA abundance declined, the linear transcriptome underwent a complete reorganization marked by the activation of novel genes. In the liver, female samples clustered by feed efficiency only when the entire linear and top-expressed circular transcriptomes were considered, respectively. In MDMs, the circular transcriptomes of control and infected samples, as well as the JD+ linear transcriptome, were dominated by donor-specific signatures. In contrast, the JD− linear transcriptome reflected MAP infection, with infection-specific structuring overriding inter-individual variation.
Conclusions
In both MDM and tissue samples, circular and linear transcriptomes follow distinct and largely independent regulatory logics. While both capture inter-individual variation, circRNA expression appears more variable and may carry fewer physiological signals, especially when no clear phenotypic signature has been detected in the corresponding linear transcriptome. These findings demonstrate that circular and linear RNAs arise from complementary and nonredundant layers of gene regulation, emphasizing the importance of analyzing both in parallel.
Background
The recent advances in high-throughput sequencing technologies and bioinformatics have considerably expanded our capacity to investigate gene transcription. Beyond the production of linear transcripts, many genes—coding or non-coding—are also capable of generating circular RNA (circRNAs) transcripts. Since the initial characterization of circRNAs in RNA sequencing data [1], the number of functions attributed to certain circRNAs has grown steadily [2, 3]. However, the functional significance of the presence of thousands of circRNAs alongside their linear counterparts in a given tissue remains poorly understood [4, 5].
In this study, we focus exclusively on exonic circRNAs, which are produced through a noncanonical splicing event known as backsplicing. These circRNAs differ from traditional linear RNAs since they lack both a 3′ poly(A) tail and a 5′ cap structure, which makes them more stable and resistant to ribonuclease-mediated degradation. Identification of circRNA relies on detecting reads that contain the circular junction. These reads, which span the junction, are aligned to the reference genome as two segments mapped in reverse order [1]. In backsplicing, a downstream splice-donor site is covalently linked to an upstream splice-acceptor site, with both splice sites coinciding precisely at known exonic boundaries [3, 6]. Specifically, the sequences of reads spanning the circular junction include the 3′ end of one exon and the 5′ start of an upstream exon [5].
Exonic circRNAs are produced by coding or noncoding poly-exonic genes [7]. Research indicates that the functions and properties of certain exonic circRNAs are not solely dependent on their sequence [8, 9, 10–11]. Exonic circRNAs are typically generated from parental genes (PGs) able to produce canonical mRNA isoforms [12]. This dual expression—producing both circular and linear transcripts—represents a significant feature of these genes. Some cases have reported circRNAs being much more abundant than their linear counterparts [1]. However, the circRNA/linear RNA ratio generally favors linear RNA [13].
The production of exonic circRNAs is largely dependent on exon/intron structure and intron content [14, 15, 16, 17, 18–19]. However, the ratio of linear to circular RNAs can be influenced by specific cellular conditions, which is supported by the findings of Liang et al. [20], who demonstrated that inhibiting or slowing down the canonical pre-mRNA processing pathway shifts the balance of protein-coding genes in favor of circular RNAs. Furthermore, it is plausible that the balance between backsplicing and canonical splicing could be regulated by changes in the physiological conditions [21, 22]. Recent studies have also highlighted that a variety of RNA-binding proteins (RBPs) can modulate the efficiency of backsplicing reactions [2, 23, 24, 25–26]. These findings suggest a comprehensive investigation into the balance between linear and circular transcripts would be a valuable research endeavor [27]. A recent study reported by Ayyildiz et al. [28] proposed that imbalanced expression of linear and circular RNAs could impact neural fitness, potentially contributing to the pathogenesis of Huntington’s disease. A preliminary comparison of circRNA expression levels assessed by total RNA sequencing (total-RNAseq) and their linear counterparts by mRNAseq—considering all genes as parents of exonic circRNAs—was conducted on a single sample and indicated that the two production pathways operate independently [15]. A novel sequencing technique, called Lexo-circSeq, is now available for analyzing the production of circular and linear transcripts from approximately 100 genes [29]. Total-RNAseq can be used to analyze both types of transcripts. As simultaneously quantifying circRNAs and linear transcripts becomes a major goal [4, 5, 22], several tools have been developed to perform this dual quantification while considering only the linear transcripts produced by circRNA PG [13, 30, 31].
The aim of this study was to investigate the behavior of circular and linear transcriptomes across different tissues and experimental conditions. We first focused on transcriptomes derived from bovine monocyte-derived macrophages (MDMs), differentiated from peripheral blood mononuclear cells (PBMCs), to examine their responses to infection by Mycobacterium avium ssp. paratuberculosis (MAP), the causative agent of Johne’s disease (JD). Using six groups of MDM samples chosen among the complete experiment implemented by Ariel et al.[32], we designed analyses based on four pairs of related sample groups. Initial efforts centered on characterizing and comparing the circular transcriptomes of the six MDM groups, while the analyses of the linear transcriptomes in MDMs were complementary, as Ariel et al. [32] had already examined these aspects in detail. To expand the scope of our investigation, we analyzed additional bovine datasets from three tissue pairs: liver (male vs. female), muscle (male vs. female), and testis (neonatal vs. post-pubertal). For each tissue, we compared the circular and linear transcriptomes using circRNAs data previously obtained [33]. These tissue analyses complemented our findings in MDMs, offering a broader perspective on circRNA and linear RNA expression across diverse biological contexts.
Results
Parallel analysis of circular and linear transcriptome of MDMs
We used 36 bovine MDM transcriptome samples originally generated by Ariel et al. [32] and Marete et al. [34], with sequencing data available at NCBI BioProject [35] under accession number PRJNA384849. These MDMs were derived from PBMCs collected from 12 cows, selected based on the presence or absence of antibodies against MAP. Animals were classified as Johne’s disease positive (JD+, cows A–F) or negative (JD−, cows G–L). PBMCs from each animal were differentiated into MDMs and divided into two series: one exposed ex vivo to MAP infection for 1, 4, 8, or 24 h and one kept as uninfected controls sampled at the same time points [32]. For the present study, we focused on MDMs infected for 24 h (Infected_24h) and uninfected controls collected at 4 h (Control_4h) or 24 h (Control_24h). This selection yielded six experimental groups, each containing six samples: Control_4h-JD−, Control_24h-JD−, Infected_24h-JD−, Control_4h-JD+, Control_24h-JD+, and Infected_24h-JD+ (see Additional file 1: ST1). For example, the three samples deriving from A-donors are A-C_4h, A-C_24h, and A-I_24h.
Circular and linear transcriptomes characterization
For the 36 bovine MDM transcriptome samples retained from PRJNA384849, an average of 56,008,965 ± 5,849,853 (mean ± SD) reads were unambiguously mapped to the cattle reference genome, out of a total of 66,567,933 ± 6,734,244 reads, yielding an average mappability of 84.18%. As no sample showed low mapping quality, all were retained for downstream analyses.
To analyze the circular and linear transcriptomes of MDMs, we used CIRCexplorer3 (CE3) [36], a combination of CIRCexplorer2 (CE2)—which identifies circRNAs originating from backsplicing—and the CLEAR module [13], which evaluates the expression of circRNAs and their linear counterparts produced by the same PG. Both expression values, FPBcirc and FPBlinear, are expressed in the same unit, FPB (fragments per billion of nucleotides), and they can be directly compared [22].
We ultimately identified 3371 exonic circRNAs in MDMs based on the analysis of approximately 2396 billion reads. Among these, 3288 had previously been reported as reliable and 13 as unreliable in the preliminary bovine circRNA catalog [37, 38]. These 13 circRNAs were removed during a final filtering step, resulting in a curated set of 3358 circRNAs derived from 1815 parental genes, consistently detected across the 36 MDM samples (Additional file 1: ST2). Notably, only 70 of the detected circRNAs had not been previously reported in bovine tissues [38].
Regarding the three expression parameters provided by CE3 (FPBcirc, FPBlinear, and CS) evaluated for 3358 circRNAs, we kept 105,612 individual results provided by the CLEAR analysis of 36 samples (Additional file 1: ST3). Since these parameters concern not only the circRNA but also its linear counterpart, it seems more appropriate to speak of Transcription_units (T_units). To obtain a complete matrix of FPBcirc, undetected circRNAs were assigned a value of zero, and the resulting cor_FPBcirc values for each of the 3358 circRNAs across the 36 samples are reported (Additional file 1: ST4). Across all 3358 circRNAs, the mean cor_FPBcirc was 0.57 FPB (SD = 0.33). Notably, the vast majority (2983 out of 3358, corresponding to approximately 89%) showed low expression levels (FPBcirc < 1).
Among the 1815 PGs analyzed, 1074 were classified as single-producer and 741 as multi-producer genes. Because all T_units from the same PG share an identical FPBlinear value, this distinction is important for accurate interpretation. The average FPBlinear across all distinct PGs was 20.11 FPB (based on data reported in Additional file 1: ST2).
Size of circular transcriptome of MDMs
When analyzing the 36 MDM samples across the 3358 bovine T_units, a total of 120,888 circRNA detection opportunities were expected (36 samples × 3358 circRNAs). However, only 105,613 individual presence calls were recorded (Additional file 1: ST3). Notably, five JD+ samples (A-C_4h, A-C_24h, A-I_24h, B-I_24h, and C-I_24h) exhibited more than 20% missing values (i.e., cor_FPBcirc = 0 in Additional file 1: ST4), whereas only one JD− sample (I-I_24h) showed a comparable pattern. This discrepancy may reflect differences in circRNA biogenesis or detection efficiency under varying infection conditions or donor backgrounds. For each MDM group derived from six donor animals, we examined the 3358 cor_FPBcirc values. On average, the coefficient of variation (CV) across these values was high (Fig. 1A), indicating substantial inter-sample variability in circRNA expression levels—particularly in Infected samples.
[See PDF for image]
Fig. 1
Circular transcriptome characterization across six MDM sample groups. A Overview of the circular transcriptome within the six MDM sample groups, including the three treatment groups from JD+ donors and the three from JD− donors. For each circRNA (whether detected or not in a sample), the mean and standard deviation (SD) of cor_FPBcirc values (reported in Additional file 1: ST4) were calculated within each group of six samples. A coefficient of variation (CV = SD/mean cor_FPBcirc) was computed for each T_unit to assess the consistency of cor_FPBcirc expression across samples. B Size of the circular transcriptome across the 36 MDMs samples
To compare the relative size of linear and circular transcriptomes within the same sample, we estimated the total number of RNA molecules for each type. The size of the circular transcriptome can be approximated by multiplying the average FPBcirc (0.57) by the number of detected circRNAs (3358), yielding ~1914 FPB. In contrast, based on data from Ariel et al. [32], who reported 9691 expressed genes in MDMs, and with an average FPBlinear of 20.11 determined above, the linear transcriptome size is estimated at ~194,886 FPB. This suggests that, overall, the circular transcriptome is approximately 100 times smaller than its linear counterpart.
Among the control samples, the size of the circular transcriptome appeared similar between JD+ and JD− groups. However, in infected samples, a notable decrease in circular transcriptome size was observed (Fig. 1A).
To explore this further, we calculated the circular transcriptome size for each MDM sample as the sum of cor_FPBcirc values and compared them across experimental conditions (Fig. 1B). In JD− donors, a consistent trend was observed: Control_4h samples exhibited the largest circular transcriptome, followed by Control_24h (except for the K-donor), and then Infected_24h samples. This pattern suggests that both prolonged culture and MAP infection are associated with a progressive decline in circular transcriptome intensity in JD− MDMs. Conversely, MDMs from JD+ donors showed more variable responses under prolonged Control conditions. The impact of MAP infection also differed: samples from C-donor, and particularly B-donor, deviated markedly from the general trend, suggesting donor-specific or atypical responses to infection that influence circRNA expression dynamics.
Comparison of circular and linear transcriptomes
To compare circular and linear transcriptomes across the 36 bovine MDM samples, we selected a set of 647 T_units, where the circRNA is consistently detected in all samples. This yielded two datasets: 647 × 36 FPBcirc values (for individual circRNAs) and 508 × 36 FPBlinear values (for the 508 linear transcripts corresponding to the 508 distinct PGs of these circRNAs). The 647 circRNAs represent a highly expressed subset, with a mean FPBcirc of 1.55—nearly three times the average observed across all 3358 circRNAs. In contrast, the 508 parental genes form a quasi-random subset, as their mean FPBlinear (22.97 FPB) is close to the overall average for the full set of 1815 parental genes (20.11 FPB) and much lower than that of the 508 most highly expressed genes (56.49 FPB).
Two hierarchical clustering analyses (HCAs) were conducted: one using the FPBcirc matrix and the other the FPBlinear matrix. The HCA based on FPBcirc values revealed two main clusters (Fig. 2A). One cluster contained only the three samples from JD+ donor A, while the second split into two subgroups—one grouping all JD− samples (donors G to L) and the other comprising the remaining JD+ samples (excluding donor A). The HCA based on FPBlinear values also identified two primary clusters (Fig. 2B). As in the FPBcirc analysis, one cluster included the three samples from JD+ donor A. In the second cluster, a subgroup of seven infected samples emerged, including six JD− donor samples and one from JD+ donor C.
[See PDF for image]
Fig. 2
Analyses of the 36 MDM samples by HCAs performed on 647 T_units, corresponding to 647 circRNAs originating from 508 distinct parental genes. Clustering was performed using the Ward method with Pearson correlation as the dissimilarity metric. A Analyses were based on FPBcirc values (647 × 36 matrix). B Analyses were based on FPBlinear values (508 × 36 matrix). Sample labels were color-coded by condition: green (Control_4h), blue (Control_24h), and red (Infected_24h). Donors JD+ and JD− were labeled A–F and G–L, respectively. Label readability has been enhanced for clarity; original plots are available in Additional file 4
According to the results obtained by HCAs, we performed the principal component analysis (PCA) separately for JD+ and JD− samples. An initial PCA was performed using 647 FPBcirc values from 18 JD− donor samples (Fig. 3A). The first two PCA dimensions capture 36.26% of the total variance, but the 18 samples remain distinguishable. Notably, samples from the same donor tend to cluster together. A second PCA with 647 FPBcirc values from 18 JD+ samples (Fig. 3C) reveals that samples from donors A and C account for over half of the total variance (Fig. 3C1). In the second and third dimensions (Fig. 3C2), samples from the same donor cluster similarly, suggesting that donor origin has a stronger influence on sample grouping than infection status, as seen with JD− samples. These PCAs, alongside the HCAs of the 36 samples, indicate that samples from donor A, and to a lesser degree donor C, diverge more from other JD+ donors in their circular transcriptome profiles. Overall, JD− samples exhibit greater homogeneity than JD+ samples in terms of circular transcriptome consistency. We performed a last PCA on circular transcriptomes using the complete matrix including 3358 cor_FPBcirc recorded for 18 JD− samples. The resulting diagram of the individual factor map (Fig. 3B) was similar to that obtained with 647 FPBcirc (Fig. 3A).
[See PDF for image]
Fig. 3
Principal component analyses (PCAs) performed on two sets of 18 samples. A Eighteen samples JD−, 647 FPBcirc. B Eighteen samples JD−, 3358 FPBcirc. C Eighteen samples JD+, 647 FPBcirc. D Eighteen samples JD+, 508 FPBlinear. E Eighteen samples JD−, 508 FPBlinear. The diagrams labeled A, B, and E pertain to principal component analyses (PCAs) conducted on 18 samples derived from donors JD−. PCAs on JD+ samples were documented in diagrams C and E. With the exception of C2, which illustrates the second and third dimensions, all other diagrams represent the initial two dimensions of the individual factor map. Diagrams D and E are relative to PCAs conducted with FPBlinear, whereas the others are relative to PCAs conducted with FPBcirc. For these PCAs (D and E), only one value per parental gene was retained to avoid redundancy. The readability of the labels on these plots has been manually improved (originals are available in Additional file 4). Control_4h, Control_24h, and Infected_24h samples were labeled in green, blue, and red, respectively. Sample positions were marked as follows: green discs for Control_4h, blue triangles for Control_24h, and red stars for Infected_24h
To investigate linear transcriptomes, two PCAs were performed on the 508 FPBlinear values: one for the 18 JD+ samples (Fig. 3D) and another for the 18 JD− samples (Fig. 3E). The PCA of JD− samples was particularly informative, with the first principal component accounting for over 71% of the total variance and clearly separating the six infected samples from the 12 controls. In contrast, the PCA of JD+ samples revealed that the linear transcriptome was more strongly influenced by donor-specific factors than by the transcriptional response to MAP infection.
Overall, multivariate analyses (HCAs and PCAs) revealed that circRNA profiles were primarily shaped by donor-specific factors—especially in MDM_JD+ —and were only marginally affected by MAP infection.
Comparison of circular and linear transcriptomes in three bovine tissues
To complement the analyses conducted on MDMs—samples derived from an extended ex vivo maturation protocol—we broadened our investigation to include bovine tissues from healthy animals. We focused on three datasets in which the circular transcriptome had already been characterized [33], reusing the existing circRNA quantifications. To this, we added an independent evaluation of linear transcript expression, enabling a parallel analysis of both transcriptomic layers within each dataset.
Specifically, we included total-RNAseq samples from the liver (six females and six males), muscle (six females and six males), and testis (T1–3 from neonatal animals and T4–6 from 13-month-old animals). The liver and muscle samples were derived from a cohort of 48 Charolais × Holstein F2 animals (NCBI BioProject [35] PRJEB34570), originally classified according to their feed efficiency [39]. For the 2021 circRNA study [33], two balanced groups had been selected from this cohort: six 18-month-old males (animals: nos. 1–3 with low feed efficiency and nos. 4–6 with high efficiency) and six 3.5-year-old females (animals: nos. 7–9 with low efficiency and nos. 10–12 with high efficiency) [39]. The testis samples (NCBI_bioproject [35] PRJNA471564; [40]) were also included in the 2021 analysis [33] and originated from six Angus bulls: three sampled during the neonatal period and three sampled at 13 months of age, shortly after the onset of puberty [41].
Distinct dynamics of the testis circular and linear transcriptomes
We began by analyzing the testis dataset, which included three neonatal and three post-pubertal samples, to examine transcriptomic changes during testis maturation. Revisiting previously characterized circRNA data [39], we found that the circular transcriptome contracted markedly with development: neonatal testes expressed 1311 highly abundant circRNAs (> 0.2 RPM on average), totaling 858 RPM, whereas post-pubertal testes expressed only 561, totaling 293 RPM (all individual data are reported in Additional file 2: ST6). Despite this reduction, these circRNAs still accounted for a substantial portion of the circular transcriptome—68% in neonates and 50% in post-pubertal samples. In parallel, the linear transcriptome expanded: the number of highly expressed genes (> 20 TPM) increased from 401 in neonatal testes to 1971 in post-pubertal ones, while overall expression levels rose more modestly, from 941,000 to 974,000 TPM (all individual data are reported in Additional file 2: ST7). This broader distribution of transcript abundance likely reflects the emergence of cell-type-specific gene expression during spermiogenesis. For example, two genes related to spermiogenesis—ENSBTAG00000019003 (TNP1) and ENSBTAG00000021493 (PRM1)—were found to be exclusively expressed in post-pubertal testes. Their expression was ranked 9th and 11th, respectively, which is consistent with their roles in late germ cell differentiation [42].
To further explore these patterns, we examined the 30 most abundant circRNAs and mRNAs in each developmental stage. The top circRNAs showed minimal overlap between stages, resulting in a combined list of 50 distinct circRNAs (Fig. 4A), consistent with two stage-specific circular transcriptomes. In contrast, the top mRNAs were more conserved, with 41 genes identified across both stages. Notably, the eight most highly expressed mRNAs were identical in neonatal and post-pubertal testes. However, many of these 41 genes displayed substantial shifts in expression levels, particularly among the 13 genes whose values ranged from 1500 to 10,000 TPM in at least one condition (Fig. 4B). Interestingly, there was no overlap between the parental genes of the top circRNAs (49 genes) and the top mRNAs (41 genes), highlighting the distinct origins of the two transcript types. For instance, the parental gene ENSBTAG00000011396 (HNRNPC) produces a circRNA that ranks 8th in expression among neonatal testis circRNAs, yet its linear transcript ranks only 485th in the same samples.
[See PDF for image]
Fig. 4
Comparison of linear and circular transcriptomes in the testis dataset. A Expression profiles of the 50 most abundant circRNAs. The 30 most highly expressed circRNAs (on average, n = 3) were identified separately in neonatal and post-pubertal testes, yielding a combined list of 50 distinct circRNAs. Each circRNA is represented by a unique color (ranked no. 1 to no. 50 by average expression in neonatal testis). Expression levels are shown for neonatal samples (left) and post-pubertal samples (right). Of the top 10 circRNAs expressed during the neonatal stage, only two (circRNAs nos. 2 and 6) are among the top 10 expressed during post-pubertal age. The individual values of the expression of these 50 distinct circRNAs were reported in Additional file 2: ST10. B Expression comparison of selected linear transcripts. The 30 most expressed mRNAs (on average, n = 3) were similarly selected in each developmental stage, resulting in 42 unique genes. Part B shows the expression of 13 genes with expression levels ranging from 1500 to 10,000 TPM in at least one stage, with neonatal samples on the left and post-pubertal samples on the right. The three largest shifts involve ENSBTAG00000042414 (SNORA23, purple), expressed only in neonatal testis, and ENSBTAG00000019003 (TNP1, light green) and ENSBTAG00000021493 (PRM1, pink), both expressed exclusively in post-pubertal testis. The individual values of the expression of these 13 genes were reported in Additional file 2: ST11. C Hierarchical clustering analyses (HCAs) of testicular samples. For each HCA, a separation into two groups of three samples was expected: neonatal testes (T1–T3) versus post-pubertal testes (T4–T6). C1 HCA based on linear transcriptomes using all genes with average expression > 2 TPM. C2 HCA based on circular transcriptomes using all circRNAs with average expression > 0.02 RPM
HCAs of the six samples using all expressed linear genes (Fig. 4C1) and all expressed circRNAs (Fig. 4C2) each successfully separated neonatal and post-pubertal groups. In the circular-based HCA, we noted that inter-individual variability was higher among post-pubertal samples than among neonatal samples. This may be related to biological expectations, as puberty is a complex and asynchronous process, making true biological replication challenging despite similar chronological ages.
Distinct clustering patterns of liver and muscle transcriptomes
We knew already that the muscle circular transcriptome is poorer (number of circRNAs, expression) than that of the liver [33] (all individual data are reported in Additional file 2: ST8 and ST9). In contrast, in both liver and muscle tissues, between 8000 and 9000 genes were found to be expressed at levels exceeding 2 TPM (all individual data are reported in Additional file 2: ST7).
When HCAs were performed using all expressed genes (> 2 TPM) or all detected circRNAs, no clear clustering patterns emerged that aligned with expected biological variables such as sex or feed efficiency. To explore possible latent structure, HCAs were then conducted separately within each tissue and sex subgroup—female liver, male liver, male muscle, and female muscle—using predefined phenotypic groupings (males: nos. 1–3 vs. nos. 4–6, females: nos. 7–9 vs. nos. 10–12) as a reference. Based on the definition of the feed efficiency phenotype established for this cohort and the animals’ documented physiological backgrounds [39], clearer clustering patterns were expected—particularly in female liver samples, where the feed efficiency phenotype was reported to be nearly detectable at the liver level at slaughter (as noted during sampling and analysis by CK).
Among the HCAs based on linear transcriptomes (Fig. 5A1, A2, A3, A4, A5), only the female liver samples (Fig. 5A4) showed the expected clustering pattern, with animals L7, L8, and L9 grouped together and distinct from L10, L11, and L12. This result was consistently observed across different analytical configurations, including those that used all annotated genes without an expression threshold (Additional file 3: page 6) and those that included all protein-coding genes regardless of expression level (Additional file 3: page 7). In contrast, reducing the number of genes considered—regardless of the selection method used—did not yield meaningful clustering results (Additional file 3: pages 2–5).
[See PDF for image]
Fig. 5
HCAs of linear and circular transcriptomes in the muscle and liver. A displays HCAs based on linear transcriptomes using all expressed genes in each tissue, i.e., all genes with expression > 2 TMP on average. B shows HCAs based on circular transcriptomes using all expressed circRNAs in each tissue, i.e., all circRNAs with expression > 0.02 RPM on average. C shows HCAs based on circular transcriptomes using a specific tissue panel. Four-hundred seventy-six circRNAs were retained for the liver (top 250 in males + top 250 in females) and 314 circRNAs for the muscle (top 150 in males + top 150 in females). For each HCA, clustering into two groups was expected. When the resulting clusters matched the expected sample’ groupings based on phenotypic information (e.g., L7–L9 vs. L10–L12 for liver females), the label background is shaded green. When the HCA led to a clustering 2 × 3, but did not match the expected groupings, the background is shaded orange-brown
In contrast, HCAs based on circular transcriptomes (Fig. 5B1B4) revealed potential 2 × 3 substructures in male muscle (Fig. 5B1) and female muscle (Fig. 5B2), though these patterns did not align with the expected feed-efficiency groupings. As improved clustering had previously been observed when focusing on the more highly expressed circRNAs [37], we performed HCAs using tissue-specific panels of the top-expressed circRNAs: 476 circRNAs were retained for liver (top 250 in males + top250 in females) and 314 circRNAs for muscle (top 150 in males + top 150 in females). Among these refined analyses (Fig. 5C), only the female liver samples displayed the expected clustering pattern (L7, L8, L9 vs. L10, L11, L12, see Fig. 5C4), while the male muscle samples again showed clustering in 2 × 3 subgroups (Fig. 5C1) that did not correspond to the expected feed-efficiency group.
Finally, only in female liver samples—and only when considering the entire linear and top-expressed circular transcriptomes, respectively—did we observe clustering according to feed efficiency. These clustering patterns appeared robust, as similar groupings were obtained using alternative dissimilarity metrics (Fig. 5A1, A2, A3, A4 vs. Additional file 3: Pages 1-8 for linear transcriptomes; Fig. 5C1, C2, C3, C4 vs. Additional file 3: Page 9 for circular transcriptomes).
Focus on impacts of the MAP infection in JD− and JD+ MDMs
To quantify the differential impact of MAP infection on the circular transcriptome in JD− and JD+ MDMs, we performed differential expression analyses (DEA) using DESeq2 on two key comparisons: Control_24h-JD+ vs. Infected_24h-JD+ and Control_24h-JD− vs. Infected_24h-JD−. To serve as negative controls, we also analyzed time-dependent changes in the absence of MAP infection: Control_4h vs. Control_24h in both JD+ and JD− groups. In all cases, a paired model was used to account for donor identity, as samples within each comparison originated from the same individuals. None of these analyses revealed any differentially expressed circRNAs.
These negative results are consistent with results observed in multivariate analyses. Given the absence of differentially expressed circRNAs at the global level, we next examined more closely the two comparisons: Control_24h-JD− vs. Infected_24h-JD− and Control_24h-JD+ vs. Infected_24h-JD+.
Controls/Infected_24h-JD−
Within these 12 samples available (6 Control_24h-JD− samples and the 6 Infected_24h-JD− samples), we identified 1524 T_units corresponding to circRNAs consistently detected in all samples. After excluding 13 T_units with low expression (mean FPBcirc < 0.2 in both groups), we retained 1511 T_units, which were derived from 1063 PGs, for comparative analyses of FPBcirc, FPBlinear, and CS values. Across all 12 samples, no correlation was observed between the 1511 FPBcirc values and their corresponding FPBlinear values (including after log2 transformation).
The FPBcirc scatter plots (Fig. 6A1) showed that data points were tightly clustered along the diagonal, symbolizing the conservation of FPBcirc under both conditions. This observation suggests a slight decrease in circRNA expression, which is consistent with the observed reduction in circular transcriptome size (Fig. 1B). In contrast, FPBlinear scatter plots, which were constructed using the log2-transformed FPBlinear values (Fig. 6A3), showed a broader distribution, indicating more variable expression. The distribution of the 9066 CS values (Fig. 6A4, log2-transformed) closely mirrored the FPBlinear scatter plot pattern. This is not surprising, as CS variability is largely driven by fluctuations in linear transcript levels—given that the circular transcriptome is approximately 100 times smaller than the linear transcriptome, and that the log2 transformation tends to favor FPBlinear over FPBcirc (see Fig. 6A2).
[See PDF for image]
Fig. 6
Comparison of transcriptomes in MDM (without/with MAP infection). A Comparison Controls/Infected_24h-JD− performed using data relative to 1511 T_units. B Comparison Controls/Infected_24h-JD+ performed using data relative to 988 T_units. A1 and B1: FPBcirc scatter plots of all available values. B3 FPBcirc scatter plots with distinction of the donor. A2 Scatter plot of all values of log2-transformed FPBcirc. A3 and B3 Scatter plots of all log₂-transformed FPBlinear values. A3 was generated using 9066 FPBlinear values, which contains some redundancies but yield the same result as using the 6378 distinct values corresponding to 1063 PGs. B3 was generated using 988 × 6 FPBlinear values, identical to a plot based on the 729 × 6 distinct values from 729 PGs. A4 Scatter plot of all values of log2-transformed CS. A5 Scatter plots comparing log2 fold changes (log2FC) for linear and circular transcripts. The x-axis represents log2(FC_FPBlinear), and the y-axis represents log2(FC_FPBcirc). Each T_unit is represented by a gray lozenge. To facilitate visualization of major transcriptomic changes, a yellow box was added highlighting T_units with variations of less than two log2 units
To provide a clear visual representation of how circular and linear transcriptomes differ in response to MAP infection, we constructed a scatter plot inspired by Gomes-Duarte et al. [43]. For each T_unit and each sample pair (Control_24h-JD− vs. Infected_24h-JD− from the same donor), we calculated two fold changes (FC): FC_FPBcirc and FC_FPBlinear, corresponding to the expression ratio (Infected vs. Control) for circRNAs and their linear counterparts, respectively. We then plotted log2(FC_FPBcirc) against log2(FC_FPBlinear), using the mean values across the six donors for each T_unit (Fig. 6A5). The resulting distribution shows a markedly narrower spread along the vertical axis (circular fold changes) than along the horizontal axis (linear fold changes), indicating that expression changes are generally more pronounced in the linear transcriptome. Notably, only linear transcripts exhibited fold changes exceeding log2FC > 2, suggesting that any biologically meaningful response to MAP infection is limited to the linear layer of the transcriptome.
All analyses of the Control_24h-JD− and Infected_24h-JD− transcriptomic data support the hypothesis that most transcriptomic changes in this context arise from the linear transcriptome. This was previously explored by Ariel et al. [31], who identified 1460 differentially expressed (DE) linear transcripts out of 9960 expressed genes. In contrast, our current analysis did not detect significant changes in the circular transcriptome under the same conditions. To further investigate potential overlap, we asked whether genes with DE-linear transcripts also produce circRNAs. Of the 1460 DE-linear transcripts reported by Ariel et al., only 97 parental genes were publicly available (Table S4 [32]). We cross-referenced these with our dataset of 1511 high-confidence T_units (Fig. 6A) and found that 6 of the 97 genes also produced circRNAs. While this proportion (6 out of 1063 circRNA-producing genes vs. 97 out of 9960 expressed genes) appears slightly lower, the difference was not statistically significant (Pearson’s chi-squared test with Yates’ correction, p-value = 0.25). These results do not support a significant under- or over-representation of circRNA-producing genes among DE-linear genes, suggesting that approximately 10–15% of circRNA host genes may also show differential expression at the linear level.
Controls/Infected_24h-JD+
For the comparison of Control_24h vs. Infected_24h groups in JD+ donors, we analyzed 988 T_units corresponding to circRNAs consistently detected across all 12 samples (6 Control_24h-JD+ and 6 Infected_24h-JD+), originating from 729 distinct PGs. Among the 5928 FPBcirc values (988 circRNAs × 6 donor pairs), we observed a strong correlation (r = 0.85; Fig. 6B1), though slightly lower than in the JD− group (r = 0.93; Fig. 6A1). We compared the distribution of these FPB value pairs to the 988 FPB value pairs obtained for each individual donor. Scatter plots for FPBcirc are shown in Fig. 6B1 and B2; FPBlinear values appear in Fig. 6B3. These comparisons revealed marked donor-specific effects. Among the six JD+ donors, MDMs from E-donor showed minimal changes in circRNA expression in response to MAP infection (trend line closest to the diagonal, Fig. 6B2), while B-donor displayed a clear reduction in expression (trend line near the horizontal axis, Fig. 6B2). These trends are consistent with the donor-specific differences in circular transcriptome size observed earlier (Fig. 1B). In contrast, the linear transcriptomes of Control_24h and Infected_24h MDMs were more consistent across donors. Nonetheless, B-donor again stood out, showing a distinct overall decrease in expression following infection (cloud of points shifted downward relative to the others, Fig. 6B3).
Contrasting inter-individual variability in circular vs. linear transcriptomes
Four HCAs were performed in parallel to compare JD− and JD+ transcriptomes. For the two HCAs based on the circular transcriptome, 1511 and 988 FPBcirc values were considered for JD− (Fig. 7A1) and JD+ samples (Fig. 7B1), respectively. For the linear transcriptome, 1063 and 729 FPBlinear values were used for JD− and JD+ samples, respectively (Fig. 7A2 and B2). As anticipated from prior analyses, only the HCA based on the linear transcriptome of the 12 samples from JD− donors produced the expected clustering pattern, clearly separating Control_24h and Infected_24h samples (Fig. 7A2). In contrast, the three other HCAs (Fig. 7A1, B1, and B2) suggest that the inclusion of samples from donor A may disrupt the clustering. When the two samples from A-donor are excluded, the variability captured by the circular transcriptome among JD+ samples (0.38, Fig. 7B3) becomes comparable to that observed among JD− samples (0.40, Fig. 7A1). In both cases, this variability appears unrelated to infection status. Similarly, the variability captured by the linear transcriptome in JD+ samples is lower (0.25, Fig. 7B4), further highlighting the limited sensitivity of the linear transcriptome to donor-specific effects. These results suggest that circular transcriptomes are more strongly influenced by the identity of the PBMC donor than linear transcriptomes, and that the A-donor is an outlier among the set of 6 JD+ donors. Additionally, dendrogram branch lengths can be used as proxies to estimate variability when clustering follows expected patterns. In the HCA of JD− linear transcriptomes (Fig. 7A2), approximately 90% of the observed variability is explained by the Control/Infected condition, while only 10% is attributed to inter-individual differences. Taken together, these observations indicate that circular transcriptomes exhibit greater sensitivity to inter-individual variation than their linear counterparts. Traces of donor identity are still present in the linear transcriptomes of MDM_JD−, but simply they are largely masked by the impact of MAP infection.
[See PDF for image]
Fig. 7
Hierarchical clustering of MDM transcriptomes by circular and linear RNA expression. A shows HCAs of 12 JD− MDM samples using 1511 FPBcirc (A1), 1063 FPBlinear (A2), and 1063 FPBcirc values, respectively. B presents HCAs of JD+ MDM samples using B1 and B3: 988 FPBcirc and B2 and B4: 729 FPBlinear values. Twelve JD+ samples were used in B1 and B2; only 10 were retained in B3 and B4 after excluding A-donor samples. Sample labels were color-coded by condition as follows: blue (Control_24h) and red (Infected_24h). Donors JD+ and JD− were labeled A–F and G–L, respectively
The HCA shown in Fig. 7A3 was built using only one circRNA per PG, selecting the one with the highest expression. This reduction in information did not negatively affect clustering quality (Fig. 7A1 and A3), which is unsurprising given the focus on highly expressed circRNAs (mean FPBcirc for 1511 T_units = 0.98; mean FPBcirc for the 1063 selected = 1.13), previously shown to enhance clustering performance based on circRNA expression [37]. The 12 circular transcriptomes of JD− donors cluster in pairs (Control_24h and Infected_24h) according to PBMC donor identity (Fig. 7A1 and A3).
Discussion
In this study, we proposed studying circular and linear transcripts in parallel using the same Total RNAseq dataset. For the analysis of MDMs, we used CIRCexplorer3 (or CE2 + CLEAR) [13], a tool specifically designed to enable direct comparison of circular and linear RNA expression through the CIRCscore metric, calculated as the ratio of FPBcirc to FPBlinear. The CIRCscore reflects the relative expression of circRNAs using the expression of their linear counterparts as a reference. However, we approached the use of this parameter with caution. Although it was developed to assess the balance between circular and linear transcript expression from each circRNA parental gene, it was not well-suited to the context of our study. Consequently, we opted not to use this parameter and instead focused directly on FPBcirc and FPBlinear values. With CIRCexplorer3, our analyses were limited to linear transcripts originating from circRNA parental genes and only in cases where the corresponding circRNAs were consistently detected across all 36 samples. This constraint posed a significant challenge, given the high variability in circRNA expression levels. Despite this limitation, however, the tool enabled us to estimate the size of the circular transcriptome overall and compare the variability captured by circular and linear transcriptomes. For the analyses performed on tissue datasets, we employed two separate tools—one dedicated to circular transcripts and the other to linear transcripts. Although this approach was significantly more time-consuming, it proved advantageous, particularly for transcriptome-wide analyses. It allowed for comprehensive coverage of the entire linear transcriptome, which would not have been possible with a unified tool focused solely on circRNA parental genes. This strategy was especially beneficial for the analysis of transcriptomes from neonatal and post-pubertal testis, where capturing the full complexity of gene expression was essential.
To investigate both the circular and linear transcriptomes in parallel, we first focused on the dataset of MDMs generated by Ariel et al. [32]. This dataset was of particular interest because it reported numerous differentially expressed genes following MAP infection in MDMs from JD− donors, but none in MDMs from JD+ donors [32]. Such a complete absence of differential expression within the same experimental context, JD+ —outside of comparisons limited to controls—is uncommon and noteworthy. A second reason for selecting this dataset was the relatively high number of samples per condition and the substantial sequencing depth available. Nevertheless, MDMs inherently carry the molecular imprint of their PBMC donors and present limitations associated with their derivation through an ex vivo maturation process, which may introduce more inter-individual variability [44].
After performing multivariate and differential expression analyses, it is definitively shown that the circular transcriptome of MDM_JD+ or JD− retained no trace of MAP infection. In contrast, a clear response is observed in the linear transcriptome between Control_24h and Infected_24h samples from JD− samples. Ariel et al. [32] identified 1460 differentially expressed linear transcripts out of 9960 in this condition, and further showed that MAP infection had little to no impact on the linear transcriptome of JD+ MDMs.
The absence of change in the circular transcriptome of MDMs, despite marked alterations in the linear one (only in MDM_JD−), underscores the lack of coordination between these two expression layers in MDM responses to infection. It suggests that biological signals, at least in this context (MDM_JD−), may primarily affect linear transcripts. However, clustering analyses show that even transcriptomes unresponsive to MAP infection (e.g., circular transcriptomes of MDM_JD− and JD+ or the linear transcriptome of MDM_JD+) still reflect genuine inter-individual variability, shaped by donor-specific factors.
This heightened sensitivity of circRNAs to inter-individual variation may be partly due to their generally low expression levels, which can hinder the reliability of quantitative assessments. Consequently, analyses such as HCA often yield more interpretable results when limited to panels of the most highly expressed circRNAs [37]. This strategy partially improves clustering performance but only partially compensates for the challenges posed by small sample sizes and subtle phenotypic differences in liver and muscle tissues. These same difficulties also affect HCAs based on linear transcriptomes in these tissues; however, restricting the analysis to the most highly expressed linear transcripts actually worsened clustering outcomes. By contrast, in MDMs, focusing on circRNA parental genes yielded clear and expected clustering patterns, likely because this gene set represents a relatively unbiased subset of expressed genes, without selection based on expression level or functional relevance.
The nature of the phenotype is indeed a critical factor. In the testis dataset, contrasting the neonatal and post-pubertal stages, both the circular and linear transcriptomes clearly reflected the strong, biologically distinct developmental phenotype. The testis dataset, despite being limited to 2 × 3 samples, yielded particularly clear insights. The linear transcriptome underwent substantial reorganization between neonatal and post-pubertal stages, while circRNA abundance significantly declined. Post-pubertal testis samples displayed a more balanced and diversified linear transcriptome, likely reflecting the incorporation of new cell-type-specific transcriptional activities, notably those related to Sertoli cells and the spermiogenic process. These distinct changes in linear and circular transcriptomes appear to be independent of each other.
Testicular samples were the most distinct, exhibiting clear transcriptomic differences, but independent, across both the circular and linear expression layers. In contrast, the feed efficiency phenotype showed no detectable transcriptomic signature in muscle tissue (both sexes) or in male liver. However, in female liver, both linear and circular transcriptomes reflected the group structure, likely due to the stronger metabolic relevance of the phenotype in this tissue and cohort. In MDM_JD-, only the linear transcriptome included a transcriptional response to MAP infection. Notably, circular transcriptomes in MDMs demonstrated greater sensitivity to inter-individual variation than their linear counterparts. Altogether, these observations support the hypothesis that circRNA expression is actively regulated, likely through mechanisms distinct from those controlling linear RNA levels. These findings align with prior studies by Ji et al. [45] and Siede et al.[46], which also concluded that circRNA regulation is largely uncoupled from that of their linear counterparts. In our dataset, no correlation was observed between circular and linear transcript expression. While Zhou et al. [47] uniquely reported a positive correlation, several other studies [15, 48, 49, 50–51] have contradicted this, suggesting independent regulation of circRNA and linear RNA expression.
Interestingly, we noted that a convergence between circular and linear transcriptomes in MDMs was observed in the Infected_24h sample from B-donor, where both transcript levels were markedly reduced compared to those of other JD+ donors. This anomaly may be due to a decrease in MDM survival for that particular sample. While large transcriptomic differences were also observed between neonatal and post-pubertal testis samples—in both the circular and linear transcriptomes—these shifts likely stem from independent regulatory mechanisms rather than coordinated changes. The linear transcriptome was profoundly re-organized, with the number of highly expressed genes nearly doubling in the 13 months testis compared to the neonatal testis. In contrast, the circular transcriptome underwent an “inverse” transformation, with a significant reduction in circRNA diversity and abundance in the 13 months testis. In 13-month testis and Infected_24h MDMs, the priority is probably to ensure the integrity and function of the linear transcriptome. One plausible hypothesis is that circRNA biogenesis diminishes when cellular metabolism undergoes reorganization. CircRNAs are believed to accumulate when the translational machinery is saturated, potentially serving as a buffering mechanism to regulate RNA availability [20]. In dynamic tissues such as the post-pubertal testis, cellular priorities appear to shift toward reorganizing the linear transcriptome and incorporating new cell types to support tissue maturation. Supporting this idea, a study in porcine pubertal testis reported that one individual with delayed puberty exhibited a higher level of circular transcripts compared to its six counterparts [15]. This observation suggests that the balance between circular and linear transcript production may be indirectly influenced by the transcriptional demands placed on mRNA during periods of developmental change.
Overall, our results suggest that differences in circular and linear transcriptomes—whether across unrelated groups or in paired comparisons—can occur independently. While other studies (e.g., in hippocampal tissues under pathological vs. healthy conditions [43] have reported coordinated changes between linear and circular RNAs, such coordination appears to be context dependent and may, in fact, be relatively uncommon. A recent study identified PTBP1 as a key regulator of circRNA biogenesis [52]. Notably, PTBP1 knockdown in K562 cells is shown to reduce global circRNA levels and impair cell proliferation, without altering the abundance of the corresponding linear transcripts. This supports the view that the expression of circular and linear RNAs is governed by distinct regulatory mechanisms.
Conclusion
Our results indicate that differences in circRNA and linear RNA expression can arise independently. In both MDM and tissue samples, circular and linear transcriptomes exhibit distinct and largely autonomous regulatory dynamics. While both RNA populations capture inter-individual variability, circRNA expression tends to be more dispersed and may carry fewer physiologically relevant signals, particularly in contexts where the corresponding linear transcriptome lacks a clear phenotypic signature.
These findings underscore the importance of analyzing circular and linear RNAs as separate yet complementary transcriptomic layers, rather than collapsing them into aggregated metrics such as the CIRCscore. A parallel analysis of both RNA populations enables the detection of transcriptomic features that would otherwise remain undetected.
Methods
Methods used for analyzing bovine MDMs
Transcriptome data from the bioproject PRJNA384849 (available at NCBI [35]) were obtained through total-RNAseq (after rRNA depletion performed with Ribo-Zero Gold kit from Illumina), generating 100-bp paired-end reads using an Illumina HiSeq 2000 platform [32, 34]. We first performed the alignment of the transcriptomic samples with STAR (v2.7.6a) [53] by applying the specific analysis parameters proposed by Cheng et al. [54], against the reference genome (Bos_taurus.ARS-UCD1.2) available at Ensembl [55]. Then, using CIRCexplorer2 (CE2) [36, 56], we performed the identification and annotation of circRNAs using the reference annotation (Ensembl release v.110) associated with the reference genome. For each detected circRNA, the CE2 output file provided genomic coordinates, the Ensembl ID of the PG, the exons involved in exonic circRNAs, and the number of reads supporting the BackSpliced Junction (BSJ). From these 36 CE2 output files, a matrix summarizing BSJ counts across all samples was constructed (Additional file 1: ST5), representing raw circRNA expression. To obtain normalized values, each of the 36 samples was further analyzed using the CLEAR module of CIRCexplorer3 [13], which reports circRNA expression as FPBcirc and provides corresponding FPBlinear values for the linear transcripts from the same PG. A list-format file was generated from the 36 CLEAR output files, containing FPBcirc, FPBlinear, and the CIRCscore (CS) for each retained circRNA (Additional file 1: ST3). The CS value, calculated by CLEAR as the FPBcirc-to-FPBlinear ratio [13, 22, 36], reflects the relative abundance of circular to linear transcripts. Occasional cases of zero FPBlinear were observed, associating with very high CS values. To prevent this bias, all T_units with FPBlinear equal to zero were excluded from the analysis.
Only circRNAs annotated by CE2 as exonic circRNAs with genomic lengths up to 200 nucleotides were retained. We compiled a comprehensive list of exonic circRNAs by first excluding rare events (when all FPBcirc available are ≤ 0.5) and then retaining only those detected in at least 20 MDM samples (Additional file 1: ST2), as recommended to ensure robustness [5, 57, 58].
Methods used for analyzing bovine tissue datasets
Testis transcriptome data from the bioproject PRJNA471564 (available at NCBI [35]) consisted of total-RNAseq of 150-bp paired-end reads generated using a HiSeq 4000 (Illumina). Muscle and liver data from the BioProject PRJEB34570 (available at NCBI [35]) were generated in parallel on an Illumina HiSeq 2500 platform with 100-bp paired-end reads. Since these two datasets were produced separately using distinct protocols and platforms (see [39] and [40] for details), only intra-dataset analyses are possible [33]. The characterization and quantification of circRNAs in these datasets were performed previously [33] and are briefly summarized here. All sequences were aligned to the bovine reference genome assembly ARS-UCD1.2. The circRNA detection was carried out using a strategy combining CE2 [56]) and CIRI2 [59], allowing us to characterize 12,589 exonic circRNAs across the 33 bovine samples considered [33]. We kept also the quantification of circRNA expression of 2021 [33] based on the number of BSJ reads provided by CIRIquant [60], corrected by the total number of uniquely mapped reads (reads per million, RPM). For the current analyses of circular transcriptomes in the liver, muscle, and testis, we retained, for each tissue, circRNAs supported by at least five BSJ reads in a sample and with an average expression level above 0.02 RPM.
For the analysis of linear transcriptomes, sequences aligned to the bovine genome using STAR were processed with RSEM (version 1.3.0 [61]) using Ensembl gene annotation v90. RSEM provides both fragments per kilobase per million (FPKM) and transcripts per million (TPM) values; TPM values were used for downstream analyses.
Statistical analyses
Hierarchical clustering analysis (HCA) and principal component analysis (PCA) were performed on transcript expression data. HCAs were performed on the Galaxy platform proposed by Sigenae [62]. All HCAs were performed via the “ward” agglomeration method as suggested by developers [63] and using Pearson’s correlations as dissimilarity metric. Data was transformed by the normalization module available on the Galaxy platform using the “log-binary” method. PCAs were also performed on this platform, with the function PCA FactoMineR. For PCA, data was transformed using the “standard-score” method. To avoid any redundancy in clustering and dimension reduction analyses using FPBlinear values, we ensured that only one FPBlinear value per PG was used throughout all HCAs and PCAs performed.
To assess intragroup variability of the expression of circRNAs, we calculated the coefficient of variation (CV), defined as the ratio of the standard deviation (SD) to the mean of the expression measured in FPB (CV = SD/FPB).
Differential expression analysis
Using the DESeq2 R package [64], differential expression analyses (DEAs) were conducted to compare a pair of subsets of the MDM samples. The reads supporting the BSJs detected by CE2 (Additional file 1: ST5) were used as input for these analyses. We used a paired statistical model in the DEA to account for the study’s paired design, where the same MDMs from each donor were subjected to both control and infected treatments (see the different donors in Additional file 1: ST1). P-values were corrected for multiple testing (p-adj) using the false discovery rate (FDR) method proposed by Benjamini and Hochberg [65]. Only circRNAs with a p-adj value less than 0.05 were considered to be differentially expressed.
Acknowledgements
MAG is funded by a predoctoral fellowship from the Junta de Castilla and León Government (Spain) and the European Social Fund. The ERASMUS+ KA131 programme has funded the training research internship of MAG at the GenPhySE group (INRAE). We sincerely thank Aroa Suárez-Vega and Juan José Arranz for their valuable contributions to the initial drafting of this manuscript.
Authors’ contributions
Conception of the work, design of the work: AR, MAG, CK, and TF. Methodology + Formal analysis + Investigation: MAG and AR with the help of LL, TF, CK, JD, and BG. Acquisition of data (Sequencing data → raw data): MAG, PD, and TF. Software: PD. MAG and AR have drafted the manuscript and all authors have substantively revised it. All authors read and approved the final manuscript.
Funding
These studies are fully associated with the FAANG initiative.
Data availability
This study used existing datasets that were already publicly available (NCBI bioproject accession number PRJNA384849 (MDMs), PRJNA471564 (Testes), and PRJEB34570 (Liver and Muscle)). All data generated or analyzed during this study are included in this published article and its supplementary information files.
Declarations
Competing interests
The authors declare no competing interests.
Abbreviations
BackSpliced Junction
CIRCexplorer2
CIRCexplorer3
Circular RNA
Values of FPBcirc corrected
CIRCscore, i.e., the ratio between FPBcirc and FPBlinear
Coefficient of variation
Differential expression analysis
Fold change
A term used to describe the expression of a circRNA in fragments per billion bases
A term used to describe the expression of a linear transcript in fragments per billion bases
Johne’s disease
Johne’s disease positive
Hierarchical clustering analysis
Monocyte-derived macrophages
MDM sample Johne’s disease negative
Mycobacterium avium subsp. paratuberculosis
RNA sequencing of mRNAs
Peripheral blood mononuclear cells
Principal component analysis
Parental gene of circRNA
RNA-binding proteins
Reads per millions
Transcripts per millions
RNA sequencing performed after ribodepletion of ribosomal RNAs
Transcription_units
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
References
1. Salzman, J; Gawad, C; Wang, PL; Lacayo, N; Brown, PO. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS ONE; 2012; 7,
2. Yang, L; Wilusz, JE; Chen, LL. Biogenesis and regulatory roles of circular RNAs. Annu Rev Cell Dev Biol; 2022; 38, pp. 263-289.1:CAS:528:DC%2BB38XhtlOku7fF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35609906][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10119891]
3. Kristensen, LS; Andersen, MS; Stagsted, LVW; Ebbesen, KK; Hansen, TB; Kjems, J. The biogenesis, biology and characterization of circular RNAs. Nat Rev Genet; 2019; 20, pp. 675-691.1:CAS:528:DC%2BC1MXhsFCgsrbF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31395983]
4. Zhang, J; Zhao, F. Circular RNA discovery with emerging sequencing and deep learning technologies. Nat Genet; 2025; 57,
5. Robic, A; Kühn, C. Circular RNA and backsplicing: unraveling the real, the misconceptions, and the unknown. Genomics Commun; 2025; 2, e004.
6. Liu, CX; Chen, LL. Circular RNAs: characterization, cellular roles, and applications. Cell; 2022; 185,
7. Robic, A; Demars, J; Kühn, C. In-depth analysis reveals production of circular RNAs from non-coding sequences. Cells; 2020; 9,
8. Busa, VF; Leung, AKL. Thrown for a (stem) loop: how RNA structure impacts circular RNA regulation and function. Methods; 2021; 196, pp. 56-67.1:CAS:528:DC%2BB3MXmtlCmtLo%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33662561]
9. Lu, D; Chatterjee, S; Xiao, K; Riedel, I; Huang, CK; Costa, A et al. A circular RNA derived from the insulin receptor locus protects against doxorubicin-induced cardiotoxicity. Eur Heart J; 2022; 43,
10. Fuchs Wightman, F; Lukin, J; Giusti, SA; Soutschek, M; Bragado, L; Pozzi, B et al. Influence of RNA circularity on target RNA-directed MicroRNA degradation. Nucleic Acids Res; 2024; 52,
11. Liu, CX; Yang, L; Chen, LL. Dynamic conformation: marching toward circular RNA function and application. Mol Cell; 2024; 84,
12. Ashwal-Fluss, R; Meyer, M; Pamudurti, NR; Ivanov, A; Bartok, O; Hanan, M et al. CircRNA biogenesis competes with pre-mRNA splicing. Mol Cell; 2014; 56,
13. Ma, XK; Wang, MR; Liu, CX; Dong, R; Carmichael, GG; Chen, LL et al. CIRCexplorer3: a CLEAR pipeline for direct comparison of circular and linear RNA expression. Genomics Proteomics Bioinformatics; 2019; 17,
14. Ragan, C; Goodall, GJ; Shirokikh, NE; Preiss, T. Insights into the biogenesis and potential functions of exonic circular RNA. Sci Rep; 2019; 9,
15. Robic, A; Faraut, T; Djebali, S; Weikard, R; Feve, K; Maman, S et al. Analysis of pig transcriptomes suggests a global regulation mechanism enabling temporary bursts of circular RNAs. RNA Biol; 2019; 16,
16. Chen YG, Kim MV, Chen X, Batista PJ, Aoyama S, Wilusz JE, et al. Sensing self and foreign circular RNAs by intron identity. Molecular cell. 2017;67(2):228–38 e5.
17. Ivanov, A; Memczak, S; Wyler, E; Torti, F; Porath, HT; Orejuela, MR et al. Analysis of intron sequences reveals hallmarks of circular RNA biogenesis in animals. Cell Rep; 2015; 10,
18. Veno, MT; Hansen, TB; Veno, ST; Clausen, BH; Grebing, M; Finsen, B et al. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol; 2015; 16, 245. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26541409][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4635978]
19. Kramer, MC; Liang, D; Tatomer, DC; Gold, B; March, ZM; Cherry, S et al. Combinatorial control of Drosophila circular RNA expression by intronic repeats, hnRNPs, and SR proteins. Genes Dev; 2015; 29,
20. Liang D, Tatomer DC, Luo Z, Wu H, Yang L, Chen LL, et al. The output of protein-coding genes shifts to circular RNAs when the pre-mRNA processing machinery is limiting. Molecular cell. 2017;68(5):940–54 e3.
21. Fumagalli, MR; Zapperi, S; La Porta, CAM. Impact of the cross-talk between circular and messenger RNAs on cell regulation. J Theor Biol; 2018; 454, pp. 386-395.1:CAS:528:DC%2BC1cXht12lsL%2FK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29964064]
22. Ma, XK; Zhai, SN; Yang, L. Approaches and challenges in genome-wide circular RNA identification and quantification. Trends Genet; 2023; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37839990][DOI: https://dx.doi.org/10.1016/j.tig.2023.09.006]
23. Chang S, Wang Y, Wang X, Liu H, Zhang T, Zheng Y, et al. HNRNPD regulates the biogenesis of circRNAs and the ratio of mRNAs to circRNAs for a set of genes. RNA biology. 2024;21(1):1–15.
24. Cao, C; Wang, C; Dai, Q; Zou, Q; Wang, T. CRBPSA: circrna-rbp interaction sites identification using sequence structural attention model. BMC Biol; 2024; 22,
25. Alfonso-Gonzalez, C; Shi, M; Gorey, S; Holec, S; Carrasco, J; Rauer, M et al. ELAV mediates circular RNA biogenesis in neurons. Genes Dev; 2025; 39, pp. 1-17.
26. Watts, ME; Oksanen, M; Lejerkrans, S; Mastropasqua, F; Gorospe, M; Tammimies, K. Circular RNAs arising from synaptic host genes during human neuronal differentiation are modulated by SFPQ RNA-binding protein. BMC Biol; 2023; 21,
27. Levacher, C; Viennot, M; Drouet, A; Beaussire, L; Coutant, S; Thery, JC et al. Disequilibrium between BRCA1 and BRCA2 circular and messenger RNAs plays a role in breast cancer. Cancers; 2023; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37046838][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10093293][DOI: https://dx.doi.org/10.3390/cancers15072176]
28. Ayyildiz, D; Bergonzoni, G; Monziani, A; Tripathi, T; Doring, J; Kerschbamer, E et al. CAG repeat expansion in the Huntington’s disease gene shapes linear and circular RNAs biogenesis. PLoS Genet; 2023; 19,
29. Naarmann-de Vries, IS; Eschenbach, J; Schudy, S; Meder, B; Dieterich, C. Targeted analysis of circRNA Expression in Patient Samples by Lexo-circSeq. Front Mol Biosci; 2022; 9, 1:CAS:528:DC%2BB38Xhsl2ju7%2FN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35755822][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9214859]875805.
30. Morlion, A; Hulstaert, E; Vromman, M; Anckaert, J; Everaert, C; Vandesompele, J et al. CiLiQuant: quantification of RNA junction reads based on their circular or linear transcript origin. Front Bioinform; 2022; 2, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36304262][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9580843]834034.
31. Di Liddo, A; de Oliveira Freitas Machado, C; Fischer, S; Ebersberger, S; Heumuller, AW; Weigand, JE et al. A combined computational pipeline to detect circular RNAs in human cancer cells under hypoxic stress. J Mol Cell Biol; 2019; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31560396][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6884703][DOI: https://dx.doi.org/10.1093/jmcb/mjz094]
32. Ariel, O; Gendron, D; Dudemaine, PL; Gevry, N; Ibeagha-Awemu, EM; Bissonnette, N. Transcriptome profiling of bovine macrophages infected by Mycobacterium avium spp. paratuberculosis depicts foam cell and innate immune tolerance phenotypes. Front Immunol; 2019; 10, 1:CAS:528:DC%2BB3cXhsVCgtrvP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31969876]2874.
33. Robic, A; Cerutti, C; Kühn, C; Faraut, T. Comparative analysis of the circular transcriptome in muscle, liver and testis in three livestock species. Front Genet; 2021; 12, 1:CAS:528:DC%2BB3MXhvVCqu7rM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34040640][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8141914]665153.
34. Marete, A; Ariel, O; Ibeagha-Awemu, E; Bissonnette, N. Identification of long non-coding RNA isolated from naturally infected macrophages and associated with bovine Johne's disease in Canadian Holstein using a combination of neural networks and logistic regression. Front Vet Sci; 2021; 8, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33969037][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8100051]639053.
35. NCBI_bioproject. https://www.ncbi.nlm.nih.gov/bioproject/.
36. Ma, XK; Xue, W; Chen, LL; Yang, L. CircExplorer pipelines for circRNA annotation and quantification from non-polyadenylated RNA-seq datasets. Methods; 2021; 196, pp. 3-10.1:CAS:528:DC%2BB3MXlsVyisbc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33588028]
37. Robic A, Hadlich F, Costa Monteiro Moreira G, Louise Clark E, Plastow G, Charlier C, et al. Innovative construction of the first reliable catalogue of bovine circular RNAs. RNA biology. 2024;21(1):52–74.
38. Robic A, Kühn C. The first reliable catalog of bovine circular RNAs. 2024;https://doi.org/10.57745/XORQHK.
39. Nolte, W; Weikard, R; Brunner, RM; Albrecht, E; Hammon, HM; Reverter, A et al. Biological network approach for the identification of regulatory long non-coding RNAs associated with metabolic efficiency in cattle. Front Genet; 2019; 10, 1130.1:CAS:528:DC%2BB3cXhs12gt7bM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31824560][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6883949]
40. Gao, Y; Li, S; Lai, Z; Zhou, Z; Wu, F; Huang, Y et al. Analysis of long non-coding RNA and mRNA expression profiling in immature and mature bovine (Bos taurus) testes. Front Genet; 2019; 10, 646.1:CAS:528:DC%2BB3cXhsFGit7w%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31333723][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6624472]
41. McGowan, M; Holland, MK; Boe-Hansen, G. Ontology and endocrinology of the reproductive system of bulls from fetus to maturity. Animal; 2018; 12,
42. Zhao, M; Shirley, CR; Mounsey, S; Meistrich, ML. Nucleoprotein transitions during spermiogenesis in mice with transition nuclear protein Tnp1 and Tnp2 mutations. Biol Reprod; 2004; 71,
43. Gomes-Duarte, A; Veno, MT; de Wit, M; Senthilkumar, K; Broekhoven, MH; van den Herik, J et al. Expression of Circ_Satb1 is decreased in mesial temporal lobe epilepsy and regulates dendritic spine morphology. Front Mol Neurosci; 2022; 15, 1:CAS:528:DC%2BB38XhvVKgs7nK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35310884][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8927295]832133.
44. Nielsen, MC; Andersen, MN; Moller, HJ. Monocyte isolation techniques significantly impact the phenotype of both isolated monocytes and derived macrophages in vitro. Immunology; 2020; 159,
45. Ji P, Wu W, Chen S, Zheng Y, Zhou L, Zhang J, et al. Expanded expression landscape and prioritization of circular RNAs in mammals. Cell reports. 2019;26(12):3444–60 e5.
46. Siede, D; Rapti, K; Gorska, AA; Katus, HA; Altmuller, J; Boeckel, JN et al. Identification of circular RNAs with host gene-independent expression in human model systems for cardiac differentiation and disease. J Mol Cell Cardiol; 2017; 109, pp. 48-56.1:CAS:528:DC%2BC2sXhtFOlsbjN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28676412]
47. Zhou, T; Xie, X; Li, M; Shi, J; Zhou, J; Knox, K et al. Rat bodymap transcriptomes reveal unique circular RNA features across tissue types and developmental stages. RNA; 2018; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30185624][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6239174][DOI: https://dx.doi.org/10.1261/rna.067132.118]
48. Rybak-Wolf, A; Stottmeister, C; Glazar, P; Jens, M; Pino, N; Giusti, S et al. Circular rnas in the mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell; 2015; 58,
49. Dong, WW; Li, HM; Qing, XR; Huang, DH; Li, HG. Identification and characterization of human testis derived circular RNAs and their existence in seminal plasma. Sci Rep; 2016; 6, 39080.1:CAS:528:DC%2BC28XitFGjsbvL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27958373][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5153637]
50. Kristensen, LS; Okholm, TLH; Veno, MT; Kjems, J. Circular rnas are abundantly expressed and upregulated during human epidermal stem cell differentiation. RNA Biol; 2018; 15,
51. Li, X; Liu, S; Zhang, L; Issaian, A; Hill, RC; Espinosa, S et al. A unified mechanism for intron and exon definition and back-splicing. Nature; 2019; 573,
52. Wang, M; Zheng, S; Zhang, Y; Zhang, J; Lai, F; Zhou, C et al. Transcriptome analysis reveals PTBP1 as a key regulator of circRNA biogenesis. BMC Biol; 2025; 23,
53. Dobin, A; Davis, CA; Schlesinger, F; Drenkow, J; Zaleski, C; Jha, S et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics; 2013; 29,
54. Cheng, J; Metge, F; Dieterich, C. Specific identification and quantification of circular RNAs from sequencing data. Bioinformatics; 2016; 32,
55. Ensembl.https://www.ensembl.org/index.html.
56. Zhang, XO; Dong, R; Zhang, Y; Zhang, JL; Luo, Z; Zhang, J et al. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Res; 2016; 26,
57. Gruhl, F; Janich, P; Kaessmann, H; Gatfield, D. Circular RNA repertoires are associated with evolutionarily young transposable elements. Elife; 2021; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34542406][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8516420][DOI: https://dx.doi.org/10.7554/eLife.67991]
58. Xu, C; Zhang, J. Mammalian circular RNAs result largely from splicing errors. Cell Rep; 2021; 36,
59. Gao, Y; Zhang, J; Zhao, F. Circular RNA identification based on multiple seed matching. Brief Bioinform; 2018; 19,
60. Zhang, J; Chen, S; Yang, J; Zhao, F. Accurate quantification of circular RNAs identifies extensive circular isoform switching events. Nat Commun; 2020; 11,
61. Li, B; Dewey, CN. RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics; 2011; 12, 323.1:CAS:528:DC%2BC3MXhtFajtLvM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21816040][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3163565]
62. SIGENAE_Galaxy_server. http://www.sigenae.org/ provides normalization, h_clust, and PCAfactoMineR tools available at https://toolshed.g2.bx.psu.edu/repos/vmarcon.
63. BIOS4BIOL. Learning clustering.http://genoweb.toulouse.inra.fr/~formation/CATIBIOS4BIOL_stats/Learning_clustering_current.pdf.
64. Love ML, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. bioRxiv. 2014.
65. Benjamini, Y; Hochberg, Y. Controlling the false discovery rate—a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol; 1995; 57, pp. 289-300.
© The Author(s) 2025. This work is published under http://creativecommons.org/licenses/by-nc-nd/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.