Prediction and quantification of circRNAs was performed with the python-based circRNA quantification and differential expression program: CIRIquant (Zhang et al., 2020). CircRNA identification was performed with the CIRI2 program that is embedded within CIRIquant. CIRI2 predicts circRNAs based on an algorithm that detects and differentiates between back-spliced junction (BSJ) reads and non-BSJ reads to reduce the number of false positives with a high level of sensitivity (Gao et al., 2015). Quantification was performed via CIRIquant and by re-alignment of RNA-seq reads against a pseudo-circular reference. Abundance of circRNA was estimated by comparison of BSJ reads and front-spliced reads (FSJ).

As a partial verification of circRNA prediction, and to assess depletion protocols for future circRNA studies, RNA of the SCs from the NCT control temperature group for both proliferation (NCT-pro-R) and differentiation (NCT-dif-R) was subjected to an RNase depletion protocol. Briefly, RNA samples were treated with RNase R enzyme to digest linear RNA molecules while retaining lariat or circular RNA structures. Digestion was performed on 4 ug of total DNAse-treated RNA using 10 units of RNase R exoribonuclease (Lucigen, Corp) following manufacturer’s protocol (incubation reaction at 37°C for 20 min followed by 65°C for 20 min). The resulting RNAs were quantified and submitted to the University of Minnesota Genomics Center for library preparation and sequencing. The libraries were sequenced on the Illumina NovaSeq SP platform. Prediction of circRNAs in depleted samples was performed with CIRIquant utilizing the “RNase R effect correction” function. The “RNase R effect correction” function uses both the depleted and un-depleted expression to estimate the before-treatment expression levels of circRNAs detected in RNase R data (Zhang et al., 2020).

Twelve expression comparisons were made between SC lines and temperature treatments. Across all comparisons, a total of 292 instances of significant differential expression involving 140 circRNAs, were found with an average of 22 circRNAs (range: 8–35) differentially expressed in each treatment comparison. Length of these circRNAs ranged from 180 to 115,731 nt, with an average length of 31,280 nt. Results for each treatment comparison are summarized in Table 2 and Supplementary Table S3. A complete list of genes associated with each DEC can be found in Supplementary Table S4.

Six DECs were found in three treatment comparisons to be differentially spliced (splice score > |.10|). All differentially spliced circRNAs (DSCs) occurred in the differentiating SCs and none were predicted in the proliferation SC comparisons (Table 2; Supplementary Table S3). Two circRNAs (circ04173; circ08492) were predicted to be down-spliced (increase in splicing towards the linear form) and one (circ01428) was up-spliced (increase in splicing towards the circular form) in the RBC2 control vs. hot treatment comparison. Three circRNAs (circ09737; circ06985; circ04168) were predicted to be up-spliced in the NCT control vs. hot treatment comparison and two circRNAs (circ04173; circ09737) were found to be predicted to be up-spliced in the NCT control vs. cold treatment comparisons (Table 2; Supplementary Table S3). Two DSCs (circ08492 and circ06985) have parent genes associated with known muscle genes, myosin light chain 10 and tropomyosin 1, respectively (Supplementary Table S4). The latter associated with actin-binding in contractile function. Three DSCs (circ04173, circ04168, circ04173) have the same parent gene, eukaryotic translation initiation factor 4 gamma 2 (eIF4G2, ENSMGAG00000005834). This gene employs various translational mechanisms to control gene expression involved in apoptosis, cell differentiation and survival, and embryonic development (Liu et al., 2023).

A total of 3,396 miRNA target binding sites (alignment score >165) corresponding to 277 turkey miRNAs was predicted for the 136 DECs with miRanda (Supplementary Table S5). For each circRNA, the miRNA with the highest alignment score (average 177.63) is given in Supplementary Table S5. Four circRNAs (circ09737; circ09900; circ09901; circ01223) had no predicted miRNA target binding sites. The number of targets per DEC ranged from 1 to 115 with an average of 24.97/circRNA. In some cases, miRNAs were predicted to bind multiple times to a single circRNA. For example, mga-miR-103 had 54 predicted target binding sites. The highest overall alignment score (200) occurred for a target binding site in circ05242 corresponding to mga-miR-26a.

Of the 277 miRNAs with target binding sites in the DECs, 43 were differentially expressed microRNAs (DEMs) in either the SC proliferation or differentiation assays (Reed et al., 2023). Eight of these miRNAs (mga-miR-206, mga-miR-2954, mga-miR-3529, mga-miR-N54, mga-miR-N40, mga-miR-N72, mga-miR-N92, and mga-miR-N185) were differentially expressed in both proliferating and differentiating SCs and had binding sites within the differentially expressed circRNAs (DECs) observed in both of the SC developmental stages.

RNAfold was used to demonstrate the potential effect of secondary structure on miRNA/circRNA interactions. Four secondary structures were predicted and examined for circ09021 (sequence data can be found in Supplementary Table S6). Included are the minimum free energy (MFE) structure, a predicted structure that does not take temperature into account and was predicted to be the same across all three temperature parameters, and three unique predicted thermodynamic ensemble (TE) structures corresponding to 37, 38°C and 39°C (TE_37, TE_38, TE_39; Supplementary Table S6; Figure 3). In general, the higher the temperature parameter, the greater the number of predicted unpaired nucleotides (Supplementary Table S6; Figure 3). Within circ09021, 11 miRNA target binding sites were predicted (alignment score >150) using miRanda. All of the predicted target binding locations were found to be unavailable for binding in the TE_37 and TE_38 structures due to loop formations. However, the TE_39 predicted secondary structure had two target binding sites exposed for mga-miR-N4 and mga-miR-N139 demonstrating the potential effect of slight temperature change on binding interactions. To investigate how slight the temperature change needed to induce a change in miRNA binding availability, RNAfold was reran to predict the secondary structure of circ09021 at 0.1°C increment changes from 38°C to 38.5°C (Supplementary Figure S1). At 38.4°C, mga-miR-N4 was able to first bind to circ09021. Three other binding sites were made available between 38°C and 38.5°C.

In all comparisons, the number of circRNAs showing increased expression was higher in cells incubated at the higher temperature. Temperature is a key regulatory factor for alternative mRNA splicing in both plants and animals (Shiina and Shimizu, 2020). Temperature-dependent alternative mRNA splicing has been demonstrated in mammals as a response to low body temperature during hibernation (Sano et al., 2015). Other effects on RNA splicing related to cellular temperature are changes in transcription rate of RNA polymerase II and RNA structure (Capovilla et al., 2015). Modifications of secondary RNA structure can affect the positioning of the spliceosomal complex which catalyzes backsplicing for the creation of circRNAs (Shi et al., 2021). These factors might contribute to the differences seen in circRNA expression in the cultured SCs. The differential expression of circRNAs with parent genes corresponding to known muscle protein genes demonstrates this potential impact. The cellular mechanisms of temperature-dependent alternative splicing continue to be the focus of ongoing research. Another potential cause to the differences in circRNA expression is that some of the circRNAs identified could simply be artifacts, though CIRIquant incorporates a multiple seed matching to help eliminate false positives (Drula et al., 2024; Gao et al., 2018)

The potential for shifts in gene expression due to thermal stress is important in poultry biology and ultimately meat production. Birds are homotherms but young birds have poor thermal regulation and they cannot maintain their body temperature when under even relatively mild thermal stress (Henrikson et al., 2018; Meltzer, 1983; Yahav, 2000). Thermal stress especially after hatch can in turn have long-term effects on the growth and structure of muscle (Patael et al., 2019; Piestun et al., 2017). Timing and degree of thermal stress are critical. Although development of satellite cells can be thermally stimulated during specific times in embryonic development, studies highlight changes in breast muscle histomorphology, metabolism, and physiology resulting from suboptimal incubation conditions (Oviedo-Rondón et al., 2020). For example, in broilers, mild thermal manipulation during early posthatch period can promote myogenic proliferation and differentiation (Halevy, 2020), whereas heat stress during the first 13 D posthatch can increase fat and collagen deposition and impair growth of the pectoralis muscle (Patael et al., 2019). Birds in production systems face potential new challenges due to climate change. These include not only higher mean temperatures, but also the increased frequency of extreme weather events, both hot and cold (Ruuskanen, Hsu, and Nord, 2021), which may in turn increase the incidence of adverse effects on muscle growth, development, and ultimately on meat quality.

Thermal comparisons within genetic lines (NCT and RBC2) found the majority of DECs were unique to genetic background with few DECs shared between lines in both proliferating and differentiating SCs (Supplementary Table S3). Differences in the DECs shared between the NCT and RBC2 lines can be attributed to underlying genetic differences from selection for commercial traits in the NCT birds. The RBC2 line been maintained at The Ohio State University, Poultry Research Center (Wooster, OH) without conscious selection for any trait since 1966 and has significantly different growth characteristics than the modern commercial bird (Lilburn and Nestor, 1991).

There is increasing evidence for the importance of RNA interactions in muscle biology. For example, our recent examination of miRNA expression in turkey muscles SCs in this same experimental system (Reed et al., 2023), detected significant expression differences among temperature treatments and between genetic lines. Including the eight miRNAs that had target binding sites with the DECs found in this study. In both proliferating and differentiating turkey SCs, significant differential expression of miRNAs important in muscle development (e.g., mga-miR-206) was found in response to thermal challenges in both turkey lines. Studies by Harding and Velleman (2016), and by Velleman and Harding (2017), suggest that miRNAs play roles in myogenic satellite cell migration. Interactions between miRNAs and circRNAs have the potential to alter these processes (Liang et al., 2017; Ouyang et al., 2018). These studies in turkey demonstrate how changes in non-coding RNAs could significantly affect cellular proliferation and differentiation, potentially affecting Pectoralis major muscle growth and development and subsequently breast meat quality.

Predictions of circRNA-miRNA interaction are often used support the hypothesis that circRNAs function as “molecular sponges”. As demonstrated here, miRanda and related programs can identify thousands of potential miRNA target binding sites associated with circRNA sequences. However, identification of potential interactions does not translate directly to downstream function and detailed biochemical experiments are necessary to verify predicted interactions. The circRNAs identified in this study and the predicted miRNA interactions provide a framework for generating future hypotheses.

Previous examinations of circRNAs have not considered how the cellular environment likely affects the function of these RNA molecules. As single stranded RNAs, circRNAs would create complex secondary and tertiary structures as a result of RNA folding. These dynamic structures respond to factors such as pH and temperature. As we have shown, a temperature increase of as little as 0.4°C induced changes in predicted secondary structure and accessibility of miRNA target binding sites (Supplementary Figure S1). These predicted changes in binding ability due to RNA structures are likely what is happening biologically within cells suggesting that even slight changes in body temperature may increase or decrease miRNA binding or binding of circRNAs to proteins. While miRNA/circRNA interactions are important to understand modulation of gene expression, changes in secondary and tertiary structures from temperature changes should also be considered important (Somero, 2018). To the best of our knowledge, this is the first time miRNA/circRNA interactions have been viewed through the lens of secondary structure.