Decoupling transcriptome layers: the distinct and variable nature of circular RNAs

Abstract

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.

Keywords: Circular transcriptome; Exonic circRNAs; Expression; Linear transcriptome; Parental genes.

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

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

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

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

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

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 log₂-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 log₂-transformed CS. A5 Scatter plots comparing log₂ fold changes (log₂FC) for linear and circular transcripts. The x-axis represents log₂(FC_FPBlinear), and the y-axis represents log₂(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 log₂ units

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