Expression of miRNAs in turkey muscle satellite cells and differential response to thermal challenge

Thermal stress alters the transcriptome and subsequent tissue physiology of poultry; thus, it can negatively impact poultry production through reduced meat quality, egg production, and health and wellbeing. The modulation of gene expression is critical to embryonic development and cell proliferation, and growing evidence suggests the role of non-coding RNAs (RNA:RNA interaction) in response to thermal stress in animals. MicroRNAs (miRNAs) comprise a class of small regulatory RNAs that modulate gene expression through posttranscriptional interactions and regulate mRNAs, potentially altering numerous cellular processes. This study was designed to identify and characterize the differential expression of miRNAs in satellite cells (SCs) from the turkey pectoralis major muscle and predict important miRNA:mRNA interactions in these developing SCs under a thermal challenge. Small RNA sequencing was performed on RNA libraries prepared from SCs cultured from 1-week-old male Nicholas commercial turkeys (NCTs) and non-selected Randombred Control Line 2 turkeys during proliferation and differentiation at the control temperature (38°C) or under a thermal challenge (33°C or 43°C). A total of 353 miRNAs (161 known and 192 novel) were detected across the sequenced libraries. Expression analysis found fewer differentially expressed miRNAs in the SCs of NCT birds, suggesting that the miRNA response to heat stress has been altered in birds selected for their modern commercial growth traits. Differentially expressed miRNAs, including those with described roles in muscle development, were detected both among temperature treatments and between genetic lines. A prominent differential expression of miR-206 was found in proliferating turkey SCs with a significant response to thermal challenges in both lines. In differentiating SCs, isoforms of miR-1 had significant differential responses, with the expression of miR-206 being mainly affected only by cold treatment. Target gene predictions and Gene Ontology analysis suggest that the differential expression of miRNAs during thermal stress could significantly affect cellular proliferation and differentiation.


Introduction
Thermal stress can have a negative impact on poultry production.Both heat and cold stress have been shown to reduce meat quality, egg production, and wellbeing, including the overall health and quality of life in production birds (Henrikson et al., 2018;Barnes et al., 2019;Patael et al., 2019;Ouchi et al., 2021).Previous studies have shown that thermal challenges can systemically create changes in live birds at the transcriptomic and physiological levels, both in specific tissues, including the important food quality tissues and muscles (Reed et al., 2017a;Reed et al., 2017b;Al-Zghoul et al., 2019;Barnes et al., 2019;Xu et al., 2021a;Reed et al., 2022a;Reed et al., 2022b).Some of these changes are tissue-specific shifts in, for example, lipid synthesis and degradation pathways and the upregulation of protein degradation mRNAs and other stress response pathways (Barnes et al., 2019;Xu et al., 2021b).Comparisons of muscles from diverse genetic lines of poultry, specifically breast muscle (pectoralis major) stem cells (satellite cells, SCs), have shown differences in their response to thermal challenges (Wilson et al., 1975;Li et al., 2018;Xu et al., 2021a;Xu et al., 2021b).As selfrenewing mesenchymal cells, SCs enable muscle hypertrophy, maintenance, and damage repair.Avian SCs are highly active in the early post-hatch period (Halevy et al., 2000;Mozdziak et al., 2002), and their activity can be altered by environmental stimuli with potential long-lasting effects on skeletal muscle growth (Piestun et al., 2013;Loyau et al., 2014).A better understanding of such responses could allow for targeted selection in breeding for creating more resilient birds.
Growing evidence suggests the role of non-coding RNAs, such as microRNAs (miRNAs), in the regulation of muscle growth and development and their response to thermal stress in animals (Andreote et al., 2014;Fu et al., 2018;Nawab et al., 2018;Sengar et al., 2018;Lang et al., 2019;Raza et al., 2021).MicroRNAs are 18-25 nt single-stranded RNAs that are thought to function primarily in posttranscriptional gene silencing by base pairing with target mRNAs, leading to destabilization, mRNA cleavage, or translational repression (Saliminejad et al., 2019).Gene silencing mediated by miRNAs plays an important role in animal development and disease (Kloosterman and Plasterk, 2006), with tissue-specific expressions being common in vertebrate development (Wienholds et al., 2005;Ason et al., 2006).Several studies have examined miRNA involvement in the skeletal muscle of poultry (Li et al., 2011;Andreote et al., 2014;Harding and Velleman, 2016;Velleman and Harding, 2017;Jebessa et al., 2018).Studies in chicken have shown that some miRNAs that are commonly differentially expressed in human muscle disorders are also differentially expressed in chicken muscledevelopment disorders (Shu et al., 2021).Other studies on chicken breast muscles have shown that miRNAs appear to be important during muscle development and the deposition of intramuscular fat, both of which can impact the final meat quality (Fu et al., 2018;Liu et al., 2021).
Little is known about miRNA expression and function in turkeys.A previous work by our group examined the role of three miRNAs (miR-16, miR-24, and miR-128) in the expression of genes essential to satellite cell function in turkeys.The inhibition of these miRNAs differentially affected the expression of syndecan-4, glypican-1, and myogenic regulatory factors, myogenic differentiation 1 (myoD) and myogenin (MYOG) (Harding and Velleman, 2016).Two of the miRNAs (miR-24 and miR-128) also played a role in myogenic satellite cell migration (Velleman and Harding, 2017).Further investigation of the miRNA expression and their functional interactions with mRNAs is needed to create a more complete picture of muscle development in production turkeys, particularly in regards to thermal stress.A critical initial step in identifying miRNA:mRNA target interactions is through miRNA characterization and computational prediction.
We have previously observed statistically significant differences in the gene expression (mRNA) of turkey p. major muscle SCs between growth-selected and non-selected birds and in response to thermal challenges (Reed et al., 2017a;Reed et al., 2017b;Reed et al., 2022a;Reed et al., 2022b).The current study was designed to identify miRNAs expressed in p. major muscle SCs, to identify promising candidates for further investigation, to characterize the differential expression of miRNAs, and to predict important miRNA:mRNA interactions in developing turkey skeletal muscle SCs.The experimental design of this study is novel, showing that the use of muscle satellite cells allows the delineation of their contribution independently from other cell types in the muscle tissue, whereby these mechanisms can be more clearly linked to satellite cell function.We hypothesized that the expression of miRNAs in turkey muscle SCs would be significantly altered by thermal challenges and would vary in cells from commercial growthselected birds compared to non-selected birds.

Materials and methods
RNA for this study was obtained from cultured SCs previously isolated from the p. major muscles of 1-week-old male Nicholas commercial turkeys (NCTs) and Randombred Control Line 2 (RBC2, representing commercial turkeys of 1966) turkeys.RBC2 turkeys were initiated in 1966 and maintained at the Poultry Research Center of The Ohio State University, Wooster, OH, without the conscious selection of any trait and were used as an important control in studies of select lines (Nestor et al., 1969).NCTs are modern meat-type turkeys obtained from Nicholas turkeys (Aviagen Group, Lewisburg, WV).

Proliferation experiment
After 24 h in the plating medium, the cells fed the growth medium and within the treatment were replicate-cultured at an experimental temperature (33 °, 38 °, or 43 °C) for 72 h with the medium being replaced every 24 h.The control temperature of 38 °C is the approximate temperature measured in newly hatched poults (38.0 °C-38.5 °C, G. Strasburg, unpublished data), and heat and cold treatments (43 °C and 33 °C, respectively) deviate from the approximate body temperature of mature turkeys (41.5 °C).These temperatures have been shown to produce significant effects on satellite cell proliferation (Clark et al., 2016;Xu et al., 2021a).At harvest, SCs were collected using RNAzol RT (Sigma-Aldrich) and stored at −80 °C until RNA isolation.

Differentiation experiment
For differentiation, the cells were cultured as previously described by Reed et al. (2017b) and Reed et al. (2022b).The cells within the treatment were replicate-plated at 38 °C (control) or at one of the challenge temperatures (33 °C or 43 °C), and the medium was changed every 24 h for the 72 h of differentiation.The cells were harvested as mentioned previously.

RNA isolation and sequencing
Total RNA was isolated from each sample by RNAzol RT (Sigma-Aldrich) extraction, DNase-treated (TURBO DNA-free TM Kit, Ambion, Inc.), and stored at −80 °C.The initial RNA concentration and quality were assessed by spectrophotometry (NanoDrop 1000), and the samples were submitted for library preparation and sequencing at the University of Minnesota Genomics Center.Each sample was further quantified by the RiboGreen assay (Invitrogen Corp.) using the 2100 Bioanalyzer (Agilent Technologies).Due to poor cell growth during proliferation at 33 °C, the RNA quantity for the RBC2 treatment group was insufficient, and this group was excluded from further analysis.For each of the remaining treatment groups, the replicate samples were prepared for sequencing (two biological replicates per treatment group).Indexed libraries (n = 22) were constructed using the Takara Bio smRNA Library Preparation Kit and sizes selected for approximately 170 bp inserts.The libraries were multiplexed and sequenced on the NovaSeq SP platform using v1.5 chemistry (Illumina, Inc.) to produce 51-bp paired-end reads (data accessioned as part of the NCBI SRA BioProject PRJNA842679).
Perlibrary FastQC reports were aggregated into a joint report for easy browsing using MultiQC 1.13 (Ewels et al., 2016).For each sequenced library, the "forward" (R1) read was used for downstream analyses.Sequencing adapters were removed from the reads using cutadapt 4.2 (Martin, 2011).The first three bases were additionally removed during trimming to remove nonbiological bases added during library preparation.Reads were removed if their lengths were shorter than 15 nt after trimming.Trimmed reads were then cleaned of ribosomal sequences using BBDuk 39.01 (https://sourceforge.net/projects/bbmap/).The reference sequences used for ribosomal depletion were large subunit and small subunit ribosomal sequences retrieved from SILVA release 132 (Quast et al., 2013).Reads were removed if they had an exact match of at least 15 nt to one of the reference sequences from SILVA.Reads that were trimmed of adapters and depleted of ribosomal sequences were used for downstream analyses.

miRNA prediction
Cleaned reads from all libraries were combined into a single file for the prediction of novel miRNAs against the turkey genome.The turkey genome assembly (GCA_943295565.1)was prepared for mapping using Bowtie 1.3.1 (Langmead et al., 2009).Characterized mature miRNAs from chicken (Gallus gallus) were downloaded from miRBase release 22.1 to use as previously known miRNAs.Novel miRNAs were predicted in the combined sequencing libraries using miRDeep2 0.1.2(Friedländer et al., 2012), and miRNA sequences were retained if their miRDeep2 score was >0.

miRNA expression profiling
Cleaned reads from each library were separately mapped to the turkey genome assembly using Bowtie 1.3.1 (Langmead et al., 2009).The options used were "-n 0 -e 80 -l 15 -m 5 --best--strata" to recreate the same parameters that were used for miRNA discovery using miRDeep2.SAM files were converted to BAM files using SAMtools 1.14 (Danecek et al., 2021).Processing of alignment files was performed in parallel with GNU parallel version 20210822 (Tange, 2018).miRNA regions identified using miRDeep2 were converted to SAF files for expression quantification.miRNA expression was quantified using "featureCounts" version 2.0.3 (Liao et al., 2014), requiring a minimum mapping quality of 10 for a read to be counted.To identify differentially expressed miRNAs (DEMs), expression values were first normalized by the library size and multiplied by a factor of 1 × 10 6 , corresponding to counts per million (CPM) mapped miRNA reads, where the library size is the total number of reads mapped to miRNA precursors.The counts matrix from featureCounts was analyzed using the "edgeR" package (Robinson et al., 2010) in the R statistical computing environment, version 4.2.2 (R Core Team, 2022).miRNAs with a low expression were filtered by removing those that did not have at least three assigned reads in at least two libraries.Global patterns in miRNA expression were assessed with principal component analysis using the prcomp function in R. Variance partitioning analyses were conducted using the "variancePartition" package (Hoffman and Schadt, 2016) in R, estimating the contributions of the incubation temperature, genotype, and an interaction between the incubation temperature and the genotypic variance in miRNA expression.Differential expression analyses were carried out with the quasi-likelihood F-test using edgeR (Chen et al., 2016).Differences were evaluated for the fold change (log 2 FC) and were considered significant at p <0.05.BioVenn (Hulsen et al. (2008) was used to create Venn diagrams.

miRNA target prediction
Potential miRNA target genes were initially predicted using TargetScan 8.0 (McGeary et al., 2019) and miRDB (Chen and Wang, 2020) using the chicken miRNA database.The given sequence differences between turkey and chicken genomes and potential binding sites were subsequently computationally predicted using miRanda v3.3a (Enright et al., 2003).Although the trio-based turkey genome assembly (GCA_ 943295565.1)offers improved coverage and assembly quality, at the time of this study, it had not been annotated in either the Ensembl or NCBI databases.Therefore, to identify the genes associated with the predicted turkey miRNAs, we used the reference UMD-5.1 assembly to align the miRDeep2predicted consensus precursor sequences via BLAST.Gene targets in turkey were predicted by aligning the miRNA sequences against all RNA transcripts in the annotated UMD-5.1 genome build (NCBI annotation 104) with position-weighted scoring, an alignment score of >150, and | energy-kcal/mol| >7.0.Enrichment tests for target genes were performed using the PANTHER Overrepresentation Test [GO Consortium release 20150430 (Mi et al., 2013); http:// geneontology.org/].GO analysis utilized the chicken (G.gallus) reference gene list with ~66% of turkey loci (Annotation 105) having ID homologs.The differences were considered significant at p <0.05. a For each library, the total number of PFs, clusters (the number of reads passing filter in millions per lane), yield (Mbases), percentage of bases with the quality score (Q) ≥30, and mean Q score are given.Frontiers in Physiology frontiersin.org05 Reed et al. 10.3389/fphys.2023.1293264Workflow availability All scripts needed to recreate the analyses described previously are available in a GitHub repository at https://github.com/TomJKono/Turkey_MSC_miRNA.

Small RNA sequencing
The results for the sequencing of the small RNAs are summarized in Table 1.The number of PF (passing filter) clusters averaged 19.52M reads encompassing an average 1,991 Mb per library.The read quality was consistently high with an average mean Q score of 32.4.Library sizes from the proliferation experiment were very similar (Figure 1A) and were averaged higher than those from the differentiation experiment where the numbers of reads were more variable among libraries (Figure 1B).Here, the number of reads was the lowest in the 33 °C treatment replicates and the highest in the 43 °C replicates.

Identification and the expression of conserved and novel miRNAs
Clean reads obtained from all sequencing libraries in both experiments were used for the detection of expressed miRNAs in the SCs using miRDeep2.The performance of miRDeep2 in the detection of known miRNAs (those identified based on the sequence comparison of their miRNA precursors with the miRBase dataset of G. gallus) and novel miRNAs is presented in Supplementary Table S1.In this study, a total of 353 miRNAs (161 known and 192 novel) were detected.The expression of putative novel miRNAs was lower than that of the known miRNAs.
Novel miRNAs were considered high-confidence if both the putative mature and star miRNAs (miRNA corresponding to the other side of the hairpin) reported by miRDeep2 were detected in at least two independent samples, having the exact same 5′-and 3′ends and allowing no mismatches.The cutoff values for confidence are somewhat arbitrary for novel miRNA predictions, but some studies suggest that a miRDeep score >1, significant RNAfold p-value, and mature reads >10 can be used as minimum values.Based on these criteria, 118 of the 192 detected novel miRNAs (61.4%) were considered high-confidence.In addition, sequencing reads were found to map to various miRBase gga-miRNAs that were not included within the "known" category.These "known miRNAs not detected using miRDeep2" were observed due to unusual precursor structures that do not fit the assumed biogenesis model of miRDeep2.The gga-miRNAs classified as "not detected" included 14 gga-miRNAs with mapped read counts >100 among the libraries (Supplementary Date Sheet S1).Although these were not included in subsequent analyses, the expression of three of these (gga-let-7b [140,731 mapped reads], gga-let-7l-2 [19,029], and gga-mir-210a [15,620]) may warrant further investigation.Using the UMD-5.1 genome assembly, 150 known and 130 novel miRNAs were uniquely mapped to the turkey genome (Supplementary File 1).Based on the NCBI UMD-5.1 annotation (v104), 181 (51.3%) of the precursor sequences (92 and 89 of the known and novel miRNAs, respectively) occurred within annotated genes.

Thermal challenge of proliferating satellite cells
Distributions of the expressed miRNAs are summarized in Figure 2. In the observed expression of 294 miRNAs, 170 were common in all treatments in both experiments (proliferation and differentiation).The expression of 39 miRNAs was low (average of 7.5 reads/treatment) and was limited to proliferating SCs.Only four miRNAs were uniquely expressed in single-treatment groups in the proliferation experiment (one each in the NCT 33 °C, NCT 38 °C, RBC2 38 °C, and RBC2 43 °C treatments), and all of them were novel miRNAs with a low expression (average of 2.6 reads/treatment).
Variation in the expression among treatment groups was visualized by the principal component analysis (PCA).In the proliferation experiment, treatment groups clustered distinctly along the first principal component (PCA1) (Figure 1C) with replicate treatment pairs occurring as nearest neighbors within the PCA space.The greatest within-treatment separation was seen for the NCT samples in the control (38 °C) treatment along the PCA2.Variance partitioning was used to estimate the contributions of the incubation temperature, genotype, and temperature × genotype interaction and found that the incubation temperature explained a greater proportion of the variation than genetic background (Figure 1E).

Identification of differentially expressed miRNAs
Normalization by the library size resulted in a counts matrix of 271 miRNAs (151 known and 120 novel) for analyses with EdgeR (Supplementary Table S2).With heat treatment (43 °C), only a single DEM (miR-206) was identified in NCT SCs in comparison to the control temperature (38 °C).The expression of this miRNA was significantly elevated by heat treatment (log 2 FC = 3.74).In contrast, heat treatment of RBC2 SCs had a greater effect on the miRNA expression where 73 DEMs (44 known and 29 novel miRNAs) were identified in comparison to the control temperature (38 °C) (Supplementary Table S3).Of these 73 DEMs, 34 were upregulated and 39 downregulated by heat treatment.Twenty nine of the 73 DEMs had |log 2 FC| >1.0, and 12 had |log 2 FC| >2.0.The greatest upregulation was observed for miR-206 (log 2 FC = 4.41), miR-N145 (3.69), and miR-N34 (2.75).The greatest downregulation was observed for the novel miRNAs miR-N77 (log 2 FC = −4.0),miR-N23 (−3.52), miR-N54 (−3.31), and miR-N82 (−2.96).Libraries for SCs proliferating at 33 °C were only sequenced for the NCT line, and no DEMs were identified in comparison to these with the control temperature.
The response of the two turkey lines was dramatically different at the control and heat treatment temperatures.At the control temperature (38 °C), four miRNAs were found to be differentially expressed between the SCs of commercial birds (NCT) and the RBC2 line (Figure 3; Supplementary Table S3).These included the known miRNAs, miR-206 and miR-184-5p, where the expression was lower in RBC2 SCs (log 2 FC = −2.29 and −1.58, respectively).The novel DEMs, miR-N96 and miR-N173, had higher expressions of RBC2 SCs (log 2 FC = 1.34 and 2.57, respectively).

Thermal challenge of differentiating SCs
A total of 255 predicted miRNAs were observed in SCs during differentiation, with only eight being uniquely expressed in the differentiating cells (Figure 2).The unique transcripts were a mix of known (six) and novel (two) miRNAs with a low average expression (3.7 reads/treatment group).Treatment groups clustered distinctly within the first two principal components (Figure 1D) and replicate samples generally clustered together with the exception of NCT 43 °C samples, which had the largest separation within the PCA space.Variance partitioning found the incubation temperature to explain the largest proportion of the variation (Figure 1F), and the genotype had a greater variance component in differentiating SCs when compared to proliferating SCs.
Five DEMs were shared between the NCT and RBC2, 43 °C vs. 38 °C, comparisons.The expression of miR-N105, miR-N140, and miR-N183 was significantly higher in heat-treated cells than that in the controls in both lines.In contrast, miR-N157 and miR-1559-5p were significantly downregulated in heat-treated cells compared to controls.
DE analysis found no significant differences in miRNA expression between RBC2 and NCT SCs at either the control temperature (38 °C) or heat treatment (43 °C).

Differentially expressed miRNAs following cold treatment
A within-line comparison found six DEMs in the NCT SCs being incubated at 33 °C relative to controls (38 °C).These included three isoforms of miR-1 and three novel miRNAs (miR-N183, miR-N140, and miR-N30) (Figure 4; Supplementary Table S5).Each of the DEMs showed higher levels of expression at 38 °C compared to 33 °C with average log 2 FC >3.0.

miRNA target predictions
Target predictions used sequences of each of the DEMs (|log 2 FC| >1.0) to query the annotated transcript sequences in the turkey genome for potential miRNA target sites.
In the heat-treated cells (43 °C), only two of the eight DEMs had |log 2 FC| >1.0.These include miR-206 (downregulated in RBC2 relative to NCTs, −1.631) and miR-N185 (upregulated in RBC2, 1.392).Targets for miR-206 are as presented previously (Supplementary Table S6), and the top targets for miR-N185 included HSPB7 (heat-shock protein family B (small) member 7) and SCN2B (sodium voltage-gated channel beta subunit 2).Given the increased number of DEMs identified in RBC2 cells under heat treatment ( 73), the number of potential miRNA interaction sites is significantly increased compared to those in NCT cells.Of the 73 DEMs, 29 had |log 2 FC| >1.0, and target site prediction identified nearly 2,500 predicted target sites within an average 1,065 genes per miRNA.Targets for the top 10 genes per miRNA are given in Supplementary Table S7.
As seen in the proliferating cells, NCT SCs showed that fewer miRNAs were significantly affected by cold treatment.In NCT cells, all six DEMs were upregulated and had |log 2 FC| >1.0, and target predictions for these miRNAs averaged 1,591.2targets within an average of 759.8 genes (Supplementary Table S11).The highest target alignment scores for these upregulated miRNAs (miR-1, three isoforms: miR-N30, miR-N140, and miR-N183) included FSD2 (fibronectin type III and SPRY domain containing 2), CSRNP1 (cysteine and serine rich nuclear protein 1), TMEM132A (transmembrane protein 132A), and the uncharacterized LOC104914368.GO analysis implicates the regulation of blood circulation (GO:1903522; 2.72×; p = 8.61E-05) as an overrepresented biological process.

Discussion
MicroRNAs are a class of small regulatory RNAs found in almost all animal species that play an important role in controlling the abundance of transcripts in the vertebrate transcriptome (Moran et al., 2017).These short RNA molecules predominantly recognize target sites in the 3′UTRs of mRNAs, typically leading to posttranscriptional repression as a means of modulating the gene expression (Simkin et al., 2020).Posttranscriptional downregulation by miRNAs can have a physiological stimulatory effect, as in the example of the Texel sheep breed where a sequence mutation produced an miR-1/ 206 binding site, leading to a muscle growth phenotype through the suppression of myostatin (Clop et al., 2006).The modulation of gene expression is critical for embryonic development and cell proliferation in poultry, and miRNAs have been reported to play important roles in these processes (Glazov et al., 2008;Hicks et al., 2008;Harding and Velleman, 2016;Velleman and Harding, 2017;Jebessa et al., 2018).The differential expression of miRNAs associated with growth traits (Li et al., 2011;Andreote et al., 2014;Ouyang et al., 2015) has been reported in chickens.In this study, an extensive set of miRNAs was characterized by small RNA sequencing of turkey p. major muscle SCs, identifying a total of 353 miRNAs (161 known and 192 novel).The presence of the known miRNA transcripts in the turkey SCs was consistent with the most abundant miRNAs observed in surveys of chicken skeletal muscles (Li et al., 2011;Ouyang et al., 2015;Khatri et al., 2018).
An expression unique to a limited set of tissues is indicative of highly specific miRNA interactions with a small set of target genes (Bassett et al., 2014).The tissue-specific expression of miRNAs may serve to broadly control translation in specific cells or developmental stages and perhaps modulate developmental fluctuation caused by the environment (Li et al., 2009).The expression of miRNAs is often elevated in specific tissues, and there is strong evidence for the action of specific miRNAs in muscle growth and development (Goljanek-Whysall et al., 2012).For example, miR-206 and closely related members of the miR-1 family are specifically expressed in mammalian muscles (Sempere et al., 2004;McCarthy, 2008;Townley-Tilson et al., 2010) and are required for proper morphogenesis during early embryonic development (Kim et al., 2006;O'Rourke et al., 2007;Ma et al., 2015).The expression of this miRNA family may also be muscle specific in poultry (Li et al., 2011).
The expression of miR-206 in mammals and chicken is enhanced by muscle transcription factors MyoD, MYOG, and myocyte enhancer factor-2 (Mef2) (Rao et al., 2006;Sweetman et al., 2008).However, its expression in mammals is inhibited by transforming growth factor-β (TGF-β) (Winbanks et al., 2011).In bovids, the inhibition of miR-206 and miR-1 was found to enhance SC proliferation (Dai et al., 2016).The downregulation of genes targeted by miR-206 is required for the transition of SCs in mice from proliferation to differentiation (Chen et al., 2006;Chen et al., 2010;Dey et al., 2011).The differential expression of miR-206 in turkeys in proliferating SCs and miR-1 isoforms in differentiating SCs suggests a significant response in SC development resulting from a thermal challenge.However, it is important to note that the functionality of this miRNA may be different in poultry.Associations between the miR-206 expression and general growth (Xu et al., 2013) and more defined traits such as birthweight (Jia et al., 2016), embryo myogenesis (Goljanek-Whysall et al., 2014), and muscle growth (Li et al., 2011) have been reported in chicken.However, few studies have characterized gene interactions with this miRNA.Search for target sites for gga-miR-206 in miRDB identified 675 predicted targets and 356 transcripts with conserved sites predicted using TargetScan.The comparison of these chicken targets with the 529 genes predicted to be targeted by miRanda in turkey found only 12 genes common to all three groups including ADPGK (ATP-dependent glucokinase), COL19A1 (collagen type XIX alpha 1 chain), FAM91A1 (family with sequence similarity 91 member A1), KTN1 (kinesin receptor), MEIS1 (Meis homeobox 1), NET1 (neuroepithelial cell-transforming 1), RAPGEF2 (rap guanine nucleotide exchange factor 2), RNF111 (ring finger protein 111), SMG7 (nonsense-mediated mRNA decay factor), TNPO1 (transportin-1), TRIM2 (tripartite motif-containing 2), and ZNF827 (zinc finger protein 827), which have various predicted cellular processes but without any notable ties to the SC function or muscle development.
Analysis of miRNAs in turkey SCs found a significant differential expression of known and novel miRNAs, both between genetic lines (RBC2 and NCT) and in response to a thermal challenge.The larger variance component attributed to temperature treatment is expected as a level of physiological response common between the genetic lines and unchanged by selection would be hypothesized.The greater number of DEMs observed in proliferating and differentiating SCs of the RBC2 line compared to the NCT suggests that miRNA response to heat stress has been altered in birds selected for their modern commercial growth traits.Previous RNA-seq studies of mRNA expression within an identical experimental system suggest that growth selection in turkeys has altered the developmental potential of SCs in commercial birds.In proliferating SCs, a greater number of differentially expressed mRNAs were observed from the growthselected NCT birds, and a pathway analysis indicated a shift toward early myogenesis (Reed et al., 2022a).In differentiating SCs, cold treatment produced expression changes in genes involved in the regulation of skeletal muscle tissue regeneration and sarcomere organization, whereas heat treatment increased the expression of genes regulating myoblast differentiation and survival, particularly in the NCT line (Reed et al., 2022b).
The function of miRNAs in gene regulation is defined by the gene or a group of genes that they target.Target predictions are important in attributing a functional consequence to miRNA differential expression.However, relying on predictions based on comparative datasets (human or chicken) is necessarily biased and highly sensitive to sequence variation due to the small interacting target sequences of miRNAs.Predictions based on the turkey genome and gene set offer a more reliable prediction and sequences.Target prediction algorithms suggest that many miRNAs may interact with a large group of genes, and this is supported by our target predictions.However, the degree to which prediction algorithms identify false positives is a concern (Pinzón et al., 2017;Fridrich et al., 2019).Therefore, the target and pathway predictions resulting from this study necessitate future validation studies to confirm miRNA-specific targets and their functions.Interestingly, three miRNAs (miR-16, miR-24, and miR-128) predicted in an earlier study (Harding and Velleman, 2016), for interacting with genes essential to the SC function (syndecan-4 and glypican-1), were expressed in the present study but were not included among the DEMs.
The interaction of miRNAs with gene targets is a function of the sequence match and accessibility of the target site, as mediated by the secondary structure of target mRNAs (Kertesz et al., 2007).Target sites for miRNAs are also subjected to variable rates of selection, and the sequence conservation of sites is a useful predictor of functionality (Krek et al., 2005).MicroRNAs appear to be under variable selective pressure ranging from strong selection acting on targets of some miRNAs to weak selection for other miRNAs that have many targets (Simkin et al., 2020).While some miRNAs and their targets are highly conserved (Chen and Rajewsky, 2006), others are genus-or species-specific (Kozomara and Griffiths-Jones, 2014).Comparative studies have shown that ancient miRNAs, those highly conserved among divergent taxa, are under stronger selection and are more broadly expressed (Simkin et al., 2020).
Studies have demonstrated that the thermal challenge affects the growth and subsequent structure of poultry breast muscles (Halevy et al., 2001;Piestun et al., 2017;Patael et al., 2019) with downstream effects on the meat quality.A thermal challenge has significant effects on SC proliferation, differentiation, and adipogenic potential with a differential impact on growth-selected lines of turkeys (Clark et al., 2016;Reed et al., 2017a;Reed et al., 2017b;Xu et al., 2021a;Xu et al., 2021b;Reed et al., 2022a;Reed et al., 2022b).The activation and proliferation of SCs is modulated by signaling molecules which direct myogenesis through signaling pathways.These processes are modulated by fine tuning gene expression, likely through RNA/RNA interactions, such as those involving miRNAs.In addition to miRNAs, we also characterize the expression of circular RNAs (circRNAs) in this same experimental system.CircRNAs are novel, single-stranded RNAs that are generated through the splicing of exonic/intronic sequences and are hypothesized to act as miRNA sinks (Wilusz, 2018).
The identification of non-coding RNA molecules provides further insight into the biological response to a thermal challenge and how selection for growth and increased muscle mass has altered this response.Our analyses identified a large number of genes and gene pathways potentially targeted by miRNAs in the turkey SCs available for future studies.The DEMs identified in this study of turkey SCs appear to be related to processes of muscle growth and development similar to their mammalian counterparts, and GO analysis suggests that the differential expression of miRNAs during a thermal challenge significantly affects cellular proliferation and differentiation.We caution that, to date, few studies have directly confirmed molecular miRNA/mRNA interactions in poultry (Velleman and Harding, 2017;Wu et al., 2019;Zhang et al., 2022), and most of the predicted gene interactions are currently based on the assumption that these RNA interactions in bird cells are similar to those observed in mammals (Goljanek-Whysall et al., 2012).There is, however, reason to assume that homologous miRNA:mRNA interactions do exist as target sites for miRNAs are amongst the most highly conserved motifs within mRNA 3′UTRs (Simkin et al., 2020).
This study identified miRNAs expressed in turkey muscle SCs, characterized their differential expression, and predicted important miRNA:mRNA interactions in turkey skeletal muscle SCs.Target gene predictions and Gene Ontology analysis suggest that the differential expression of miRNAs during a thermal challenge could significantly affect SC proliferation and differentiation.The distribution of DEMs suggests that selection for commercial production traits has altered the miRNA expression, providing new hypotheses for future research.

FIGURE 1
FIGURE 1 Distribution of sequencing reads (millions of fragments) in sequencing libraries of treatment groups of (A) proliferating and (B) differentiating muscle satellite cells.PCA plots of normalized read counts for (C) proliferating and (D) differentiating satellite cells.Sample-to-sample distances (within and between treatments) are illustrated for each treatment sample on the first two principal components.The samples are plotted according to the treatment.The distribution of sample variance by the treatment factor: genotype (line), temperature, interaction, and residual for (E) proliferating and (F) differentiating satellite cells.

FIGURE 2
FIGURE 2UpSet plot(Conway et al., 2017) of the expressed miRNAs.For inclusion, miRNAs must first have at least three assigned reads in at least two libraries and a treatment group average number of reads of >2.0.The horizontal bars on the left indicate the number of miRNAs expressed in each treatment.Individual points in the matrix represent miRNAs specific to each treatment, and the lines between points represent the miRNAs common to different groups.We excluded 70 single miRNAs that individually had unique group distributions.The vertical bars above indicate the number of miRNAs specific to or common to different treatments.The distribution of miRNAs by the number of treatment groups in which they were included is given above.

FIGURE 3
FIGURE 3 Distribution of DEMs during the proliferation of cultured turkey p. major SCs.For each temperature comparison, the DEMs with FDR p-value <0.05 that were shared or unique to each line (RBC2 and NCT) are indicated in the Venn diagram.The circle size is proportional to the number of DEMs.

FIGURE 4
FIGURE 4 Distribution of DEMs during the differentiation of cultured turkey p. major SCs.For each temperature comparison, the DEMs with FDR p-value <0.05 that were shared or unique to each line (RBC2 and NCT) are indicated in the Venn diagram.The circle size is proportional to the number of DEMs.

TABLE 1
Summary of RNA-seq data used for miRNA discovery and expression analyses a .