Transcriptome profiling reveals stress-responsive gene networks in cattle muscles

Animals and samples

The experiment was conducted with 32 pure Normand cull cows of 48–60 months of age purchased from different private farms in western France. The cows were not pregnant or lactating and had a medium fatness score. They were housed in the experimental farm of the INRAE research centre (UE Herbipôle-Low Mountain Ruminant Farming Systems Facility; DOI 10.15454/1.5572318050509348E12) as described by Gobert et al. (2009) and Delosière et al. (2020).

The cows (mean live weight of 642 kg) received a straw (30%) and concentrate (70%) based diet supplemented with lipids (40 g oil/kg diet DM) obtained from extruded linseeds (60%) and rapeseeds (40%) for 101 ± 3 d. For one group of cows, this diet was the control diet described in Delosière et al. (2020). For the other group, the diet was supplemented with vitamin E (155 IU/kg diet DM) and plant extracts rich in polyphenols (INRA patent #P170-B-23.495 FR; c7 g/kg diet DM, respectively; EP diet). Supplementation with vitamin E (lipophilic antioxidant) and plant extracts rich in polyphenols was used to protect against beef lipoperoxidation in a previous study on sheep (Gladine et al., 2007). The cows received a morning and evening meal representing a daily quantity of 1.8 kg of concentrate and 0.8 kg of hay. They had free access to water. The cows received an isoenergetic and isonitrogenous ration for a mean daily body weight gain of 1.6 kg over the 101 ± 3 d finishing period. Animals were housed in groups of four in 6 × 6 m pens with straw bedding, according to a balanced design relative to feeding treatments. The cows were housed in a freestall barn equipped with individual feed bunks and automatic gates. They were individually offered their appropriate allowance of concentrates and straw each day. Feed intake was calculated daily and adjusted periodically to ensure a daily gain of 1,150 g/d. The finishing period of 101 d, slightly longer than the French standards for cattle (approximately 70 d), was chosen to achieve good production conditions and to facilitate experimental organization.

Cows were finished and slaughtered under conditions of limited stress (n = 16) or physical and psychological stress (n = 16) as described in Bourguet et al. (2010). Groups included the same numbers of animals for each diet. For the limited stress condition, the cows were transported with a nonexperimental conspecific to avoid social isolation stress and were handled calmly. Specifically, for the stress condition, each cow was individually transported in a lorry (social isolation) towards an unfamiliar farm, where it was taken through a labyrinth by two purposefully noisy experimenters over a period of 30 min. The cow was then transported for 15 min to the experimental slaughterhouse. All cows were slaughtered without any electrical stimulation in compliance with INRAE ethical guidelines for animal welfare at the experimental slaughterhouse of INRAE. The cows were stunned by a captive bolt gun and exsanguinated, as performed in French commercial slaughterhouses. Carcasses were stored in a chilling room (4 °C) for approximately 45 min following exsanguination. Carcasses were sold for human consumption as in any controlled commercial slaughterhouse. Stress status was evaluated through the plasma and urinary cortisol and heart rate as described in Bourguet et al. (2010). Muscle samples from the M. semitendinosus (ST) (a himdmuscle involved in locomotion) and the M. longissimus thoracis (LT) (a support muscle for the body) were excised within 15 min after slaughter. They were immediately frozen in liquid nitrogen and stored at −80 °C until RNA extraction.

As indicated in Bourguet et al. (2010) and Delosière et al. (2020), experimental procedures and animal holding facilities respected French animal protection legislation, including licensing of experimenters. The procedures were controlled and approved by the French Veterinary Services (agreement B63 345 17). The animal experimental design was described in and registered in the research unit quality management database.

Muscle transcriptome analysis

Transcriptome analyses were carried out with Agilent oligonucleotide microarrays designed with 10,064 probes (including 1,614 control probes) for 4,210 bovine genes, including more than 3,000 specific muscular genes. The microarray was designed for monitoring transcriptional changes in genes involved in muscle growth (including energy and protein metabolism), carcass composition, fat metabolism and beef quality (including marbling). The microarray was first described in Hocquette et al. (2012a) and subsequently used in Costa et al. (2018).

RNA extraction, RNA quality checking and quantification, target amplification and labelling, hybridization with the probes, and extraction of fluorescent hybridization signals were performed as previously described in Hocquette et al. (2012a).

Data were preprocessed by Feature Extraction 10.1 software (Agilent Technologies, Santa Clara, CA, USA) for all samples and probes. The probes not meeting the quality criteria (saturation and uniformity of spots, intensity above background noise, etc.) were filtered out. Each array was normalized by dividing the raw intensity values of its probe by the median intensity of the control probes of the array. Each probe intensity was then normalized by dividing its raw value by the median of the corresponding probes from all arrays. After removing the probes with missing values, a log₂ transformation was applied to the data.

The transcriptomic data were submitted to Gene Expression Omnibus (GEO) under the accession number GSE119912. Differential analyses were conducted via linear modelling with the diet supplementation* period * stress interaction factor to explain the probe levels. An empirical Bayes method was used to moderate the standard errors of the estimated log-fold changes using the R/Limma package (http://bioinf.wehi.edu.au/limma/) as described in Smyth, Yang & Speed (2003) with a Benjamini and Hochberg multiple testing correction (Benjamini & Hochberg, 1995). The genes for which at least 80% of the probes were consistently different at the adjusted p value 10% (i.e., estimated rate of true positives in the probe list of 90%) were retained and considered differentially expressed genes (DEGs). All probe ratios were found to be consistent for each DEG, meaning that for one gene, all probe ratios were lower than 1 or greater than 1.

Gene Ontology enrichment

Functional enrichment according to Gene Ontology Biological Process (GO BP) was assessed by submitting lists of accession numbers (for DEGs) or gene names (for common regulators and targets of the DEGs) to ProteINSIDEv2 (Kaspric et al., 2015, https://umrh-bioinfo.clermont.inrae.fr/ProteINSIDE_2/). This workflow enables the analysis of lists of protein or gene identifiers from ruminant species and gathers biological information provided by functional annotations, putative protein secretion and protein interactions. It queries the g:Profiler database based on the most complete information available for Bos taurus. The list of array probes was used as a background list for enrichment analysis of the DEG lists. The GO enrichment test was declared significant at a Benjamini–Hochberg FDR < 0.08 (i.e., estimated rate of true positives in the gene list of 92%). The results are expressed as –log₁₀ (p value) on the graphs.

Identification of cis-transcriptional modules

Promoter sequences were extracted using the program Gene2Promoter (version 3.4.1; Genomatix Software Suite, Munich, Germany, www.genomatix.de) using the default settings, 500 bp upstream and 100 bp downstream of the transcription start site. We selected bovine promoters with at least one relevant transcript and preferentially a high quality level (experimentally verified 5′ transcript or with 5′ end confirmed by PromoterInspector prediction) and for whose number of conserved orthologous promoters was at least 50% of loci. This was performed for DEGs and for all of the genes in the microarray. As the coregulation of mammalian genes usually depends on a combination of TFs rather than individual TFs alone, cis-acting regulatory elements are often organized into frameworks of motifs called cis-transcriptional modules. The selected promoters were submitted to the ‘ModelInspector’ task of GEMS Launcher (version 4.1; Genomatix Software, Munich, Germany, www.genomatix.de) to search for cis-transcriptional modules. For this purpose, the promoter sequences of the genes were scanned for matches to the Promoter Module 5.4 Library (Vertebrate Module section). Fisher’s exact test was then used to identify overrepresented cis-transcriptional modules in the DEG set compared to the total gene set of the microarray.

Construction of interaction networks

Network analysis was performed with Pathway Studio software version 12.0.1.9 using Elsevier’s Resnet Mammal DataBase (Ariadne Genomics, Rockville, MD, USA). Gene interaction networks were built with the DEGs targeted by the most represented cis-transcriptional modules for each muscle (targeting at least 5 and 4 DEGs for LT and ST, respectively) and with the DEGs targeted by the nine overrepresented cis transcriptional modules common to both muscles, generating two muscle-specific networks and one common network of stress response. For each set of genes, to reconstruct the network, Pathway Studio searches known relations between the genes and adds regulators and the expression targets common to them. Filters were applied to identify only key expression regulators and targets of each network. To be added to the network, target genes must be linked to a minimum of three bibliographic references and have at least six known relations in the Pathway Studio Database. For regulators, three bibliographic references and two known relations are also needed, except for the regulators of the nine common modules for which a cut-off of five relations was chosen.

Finally, Venn diagrams were used to identify the major regulator genes and major target genes among the DEGs targeted by cis-transcriptional modules specific to LT and ST and in both muscles. Subnetworks between DEGs and their major regulator genes and between DEGs and their major target genes were extracted.

Validation of differential expression levels

A RT–qPCR assay was performed on the LT samples of 10 animals/group for 4 genes (ATF3, CEBPD, SMAD7 and FOS) with the StepOne Plus™ Real-Time PCR System using the Power SYBR1 Green master mix (both Applied Biosystems, Foster City, CA, USA). The GeNorm algorithm (Vandesompele et al., 2002) was used to determine the optimal number of reference genes required to effectively normalize the qPCR data. Four housekeeping genes were selected: UXT, MRPL39, CLN3 and TOP2B. The primer sequences (Table S1) with an annealing temperature of 60 °C were designed using Primer3 software. qPCR was performed using a StepOnePlus thermocycler (Applied Biosystems, Foster City, CA, USA). The PCR efficiency of each primer pair was tested with a 10-fold dilution series of purified cDNA. Each reaction was subjected to melting curve analysis to ensure the specificity and integrity of the PCR product. Student’s t test was used to test the significance of the difference between the limited stress and stress groups.

Identification of genes corresponding to DEGs in quantitative trait loci (QTLs)

A query of genetic information from the lists of the DEGs and the common regulators and targets of the DEGs was performed with the QTL module included in ProteINSIDEv2. Briefly, each DEG was searched on NCBI to retrieve the location of the corresponding gene on the genome. Then, this location was compared with the QTL positions in the QTL database “AnimalQTLdb”. The location on the genome must be included entirely in a QTL to consider that the DEG is mapped in the QTL.

Transcriptomic profiles

We recorded changes in gene expression profiles in the LT and the ST. Individual data are available in the GEO repository under accession number GSE119912. No effect was detected for the EP diet, the stress*EP diet or the stress*diet in either muscle, but an effect of stress was detected (P < 0.1). In the stressed cows compared to cows handled with limited stress, microarray analysis revealed changes in the abundance of 67 transcripts in the LT (including 43 up- and 24 downregulated transcripts; P < 0.1, Table S2) corresponding to 62 DEGs with unique gene names (Fig. 1). In the ST, changes in 36 transcripts were detected (including 33 up- and 3 downregulated transcripts; P < 0.1, Table S2) corresponding to 32 DEGs with unique gene names (Fig. 1). Among the differentially expressed transcripts, 27 were common to both muscles, corresponding to 24 unique gene names (Fig. 1). They included eight known transcription factors (TFs): SMAD7, ETS2, MYOG, ATF3, HES6, CEBPD, HEYL, and FOS (Table S2). In addition, muscle-specific genes were differentially expressed according to the stress status (38 in the LT and 8 in the ST; Fig. 1). These genes included four TFs (MYOD1, MYF6, CEBPB, and HES1) and one transcription cofactor (MED23) in the LT and a transcriptional activator (ZNF750) in the ST. The differential abundance of four TF transcripts (ATF3, CEBPD, SMAD7 and FOS) was checked by qPCR experiments in the LT and confirmed the observed changes, as illustrated in Table 1.

Lists of DEGs according to stress status were submitted to biological information mining through Gene Ontology (GO) term enrichment compared to the background list of the microarray (Data S1). In LT, 9 GO biological process (GO BP) terms were enriched (P < 0.08). In the ST, 26 GO BP terms were enriched (P < 0.08). As illustrated in Fig. 2, 9 of these GO terms identified in both the LT and ST: regulation of gene expression (23 genes in the LT, 13 in the ST), transcription by RNA polymerase II (16 genes in the LT, 12 in the ST), regulation of transcription by RNA polymerase II (16 genes in the LT, 12 in the ST), regulation of biosynthetic process (21 genes in the LT, 14 in the ST), regulation of cellular biosynthetic process (20 genes in the LT, 13 in the ST), regulation of macromolecule biosynthetic process (20 genes in the LT, 13 in the ST), regulation of cellular macromolecule biosynthetic process (20 genes in the LT, 13 in the ST), skeletal muscle cell differentiation (4 genes in the LT, 3 in the ST), and muscle organ development (8 genes in the LT, 5 in the ST). For each considered GO BP, the list of genes included both common and muscle-specific DEGs. There were no significant GO cellular component terms at P < 0.08 in both the LT and ST muscles (data not shown). In the LT, 18 GO molecular factor (MF) terms were enriched (P < 0.08). In the ST, 19 GO MF terms were enriched (P < 0.08). Sixteen of these GO terms were enriched in both the LT and ST. They mainly refer to DNA-binding and transcription regulator activity (Supplemental Dataset).

Cis-transcriptional modules

A promoter analysis was performed with Gene2promoter of the Genomatix Software Suite to identify common TF binding sites in the promoter regions of genes–called cis transcriptional modules–that may account for coregulation among differential transcripts. For 52 of the DEGs in the LT, 168 promoters were retrieved from the Genomatix Promoter Database, of which 111 were selected for further analysis (Table S3). ModelInspector enabled the retrieval of 288 different cis-transcriptional modules (at 1,378 locations). For 28 of the DEGs in the ST, 84 promoters were retrieved. Of these, 57 promoters were further analysed with ModelInspector, and 201 cis-transcriptional modules were found (on 675 match positions). The same analysis was performed for all the genes represented on the microarray. As illustrated in Table S3, 24 cis-transcriptional modules were identified as overrepresented DEGs compared to the genes represented on the microarray (P < 0.1) in the LT and 25 in the ST. Nine of the modules were overrepresented in both muscles. The cis-transcriptional modules and the DEGs targeted by these modules in each muscle as identified by ModelInspector are listed in Table 2. Cis-transcriptional modules with binding sites for TFs in the ETS family and SP1 family had a high occurrence in the promoters of the DEGs in both muscles.

Interaction networks and identification of regulators and main targets of DEGs

Finally, with Pathway Studio 2, we constructed interaction networks between the DEGs targeted by the overrepresented cis-transcriptional modules for each muscle and between the DEGs targeted by the 9 overrepresented cis-transcriptional modules common to both muscles. We thus generated muscle-specific networks and one core network of stress responses. Then, using the Pathway Studio 2 database, we searched for the main regulators and the main targets of the 3 networks (Data S2, sheets 1–6). We next identified the similarities among the lists obtained from these datasets to identify the key common regulators and targets (Data S2, sheets 7–8). Ten main regulators of stress-responsive genes were identified: AKT1, EGF, HIF1A, IFNG, IL1B, INS, MAPK1, MAPK14, TGFB1, and TNF. GO mining showed that they were highly related (P < 0.001) to the regulation of oxidoreductase activity and specifically the regulation of monooxygenase activity, inflammatory response, immune system, and carbohydrate metabolism (transport and metabolic processes) (Data S2, sheet 7). A list of 14 main targets of the DEGs was identified: BCL2, BGLAP, CDKN1A, COL3A1, ERBB2, FN1, ICAM1, IL6, MMP2, PPARG, SELE, SLC2A4, TLR4, and VEGFA. GO mining showed that they were highly related (P < 0.001) to response to cytokine, response to oxygen levels, response to glucocorticoids, response to stress, immune system, angiogenesis, and carbohydrate homeostasis (Data S2, sheet 8). The networks of the DEGs and the common regulators and of the DEGs and common targets are presented in Fig. 3. The list of DEGs (LT-specific, ST-specific and common DEGs), the main regulators and the main targets of the DEGs are summarized in Table 3. Twenty-four of the corresponding genes were located in QTLs linked to meat quality: shear force (CDIPT, CEBPD, DNAJB4, GPAM, RAB3IL1, MAPK1, and TLR4), muscle compression (ADRB2), tenderness score (ADRB2, CDIPT, RAB3IL1, and IFNG), muscle pH (DLL4, ERBB2), juiciness (ATP1B1, DFFB, RAB3IL1, SELE), and marbling (CEBPD, DLL4, ERBB2, GADD45, ICAM1, IL1B, IL6, LEAP2, MYF6, PDK4, PMP22, SMAD7, TNF). Some genes were also mapped in QTLs associated with marbling score (CEBPD, DLL4, DNAJB4, GADD45A, LEAP2, MYF6, PDK4, PIGM, PMP22, RAB3IL1, and SMAD7) and with lipid class contents (e.g., monounsaturated fatty acid content and conjugated linoleic acid content: ADRB2; omega-6 to omega-3 fatty acid ratio and palmitic acid content: LCAT; trans-11, cis-15-C8:2 fatty acid content: GADD45A).

Transcriptional response to stress

Understanding how preslaughter stress impacts muscle physiology would provide information for the management of beef quality, especially tenderness. In this study, we compared the muscle transcriptional profiles of cows exposed to preslaughter emotional and physical stress with those of control cows handled with limited stress. We hypothesized that this approach may be useful for investigating the molecular mechanisms of the stress response and the potential impact on meat quality. We identified changes in the abundance of several gene transcripts in two muscles of the cows exposed to stress. We found evidence of a common transcriptional response in both muscles, albeit with different metabolic types and activities, even though some muscle-specific DEGs were detected. Notably, there was a core stress response in both muscles, as shown by common DEGs and common GO terms (mainly related to the regulation of gene expression and muscle development). The highest number of DEGs was detected in the LT. This may be related to the high oxidative metabolism of the LT (Hocquette et al., 2012b), which makes it more prone to changes in oxidative status and therefore susceptible to stress. Muscle gene expression in response to stress likely also depends on sex (Oster et al., 2014), the nature and intensity or duration of the stress, and breed, which may explain some differences in the results of our study and a previous study on Angus animals (Zhao et al., 2012). Nevertheless, the regulation of genes involved in carbohydrate, lipid, and protein metabolism is likely to occur in many cases, as observed in this study, as well as in other studies on cattle (Buckham Sporer et al., 2007; Zhao et al., 2012) and pigs (Davoli et al., 2009).

Newly translated transcription factors and their related biological pathways

While the short-term response to stress may be primarily driven by changes in protein phosphorylation (e.g., reversible phosphorylation (Mato et al., 2019)), as well as enzyme activity or protein abundance, our study provided convincing evidence that the response to stress includes a transcriptional component, as previously reported in two studies (Davoli et al., 2009; Zhao et al., 2012). Indeed, functional annotation of the lists of DEGs revealed the enrichment of GO terms related to the regulation of gene expression and transcription. It is well accepted that the primary response to stress involves the activation of pre-existing TFs by phosphorylation (Sabban & Kvetňanský, 2001). Our data indicate that newly translated TFs may also relay the stress response, as 11 of the DEGs encode TFs. Eight of these TFs were identified in both muscles, of which some were detected as nodes in the molecular networks associated with the response to stress. The majority of the differentially expressed TFs were upregulated, except for two muscle regulatory factors (MYOG and MYOD1) and a transcriptional repressor (HES1). MYOG and MYOD1 are basic helix-loop-helix family TFs essential for myogenesis, including during the regenerative process (Zammit, 2017). HES1 is a downstream target of Notch (Borggrefe & Oswald, 2009). It is also a master regulator of glucocorticoid receptor-dependent gene expression. It is silenced by the primary stress hormones glucocorticoids (Revollo et al., 2013). The observed downregulation of HES1 was not surprising since stressed cows showed high plasma and urinary cortisol levels (Bourguet et al., 2010). Of the upregulated TFs, 4 were basic leucine zipper (bZip) TFs: FOS, ATF3, CEBPB, and CEBPD. ATF3, a member of the mammalian cAMP-responsive element-binding protein (CREB) family, is induced by various stresses. ATF3 is a sensor for a wide range of conditions and modulates the immune response, atherogenesis, cell cycle, apoptosis, and glucose homeostasis (Jadhav & Zhang, 2017). ATF3 has been considered an adaptive response gene with a dual mode of action to activate (as a homodimer) or repress (as a heterodimer) target gene expression. It was proposed that ATF3 functions as a “hub” of the cellular adaptive-response network that helps cells adapt to disturbances of homeostasis (Hai, Wolford & Chang, 2010). ATF3 was also found to be differentially expressed following acute stress induced by surgery in Angus beef (Zhao et al., 2012). The bZip proteins CEBPB and CEBPD are members of the C/EBP family, which participates in a number of biological responses, including energy metabolism, cell proliferation and differentiation, and immune responses (Ramji & Foka, 2002). Their binding sites are found in the regulatory regions of a large number of acute phase proteins. A dual role was proposed for the C/EBP proteins as mediators of both inflammatory responses and glucocorticoid effects (Nerlov, 2007; Roos & Nord, 2011). CEBPD expression is induced by inflammatory effectors and hypoxia and promotes proinflammatory signalling and hypoxia adaptation (Balamurugan & Sterneck, 2013). CEBPB was also recently identified as a novel regulator of satellite cell homeostasis that promotes differentiation at the expense of self-renewal (Lala-Tabbert et al., 2016).

Overrepresented TF binding sites in the promoters of DEGs and related biological pathways

Several cis-transcriptional modules were located in the promoters of the DEGs. Modules common to both muscles were detected mainly in the promoters of common DEGs, while muscle-specific cis-transcriptional modules were detected in the promoters of muscle-specific DEGs. However, some specific cis-transcriptional modules were detected in the promoters of common DEGs. FOS was targeted by 4 common modules in both muscles plus 1 specific module in the ST. ATF3 was targeted by 1 common module in both muscles and by 3 specific modules (1 in the LT and 2 in the ST). Examination of cis-transcriptional modules of DEGs from both muscles revealed that binding sites for the transcription factor SP1 and for members of the ETS family are often included in those modules. SP1 is ubiquitously expressed and, in addition to functioning as a ‘housekeeping’ TF, may be a key mediator of gene expression induced by insulin and other hormones (Solomon et al., 2008). ETS1 is a highly conserved TF that controls the expression of cytokines, chemokines and angiogenesis factors (Russell & Garrett-Sinha, 2010). ETS binding sites were found in the promoter of common differentially expressed TFs as well as in the promoter of 11 of the 15 LT-specific DEGs.

Other biological pathways related to the response to stress

Our study provided additional evidence that the response to stress interacts with the immune response, inflammatory response, and chemotaxis, as well as the production of interleukins (IL-16 in the LT; IL-1 B and IL-6 as main regulators and targets of the DEGs; and IL-10 and IL-13 in the list of common targets of the DEGs). This is consistent with previous studies examining the response to stress in livestock animals: A transcriptional shift in acquired and innate immunity pathways was reported in the peripheral blood of psychosocially stressed pigs (Oster et al., 2014). Amplified inflammatory activity was also detected in blood neutrophil expression in young bulls following truck transportation for 9 h (Buckham Sporer et al., 2007) and in the LT muscle of Angus beef exposed to acute stress induced by surgery (Zhao et al., 2012). Moreover, a conserved transcriptional response to chronic social stress involving increased expression of proinflammatory genes (including IL-6 and IL-8) has been reported in blood leukocytes (Powell et al., 2013) in mice and humans. In our study, changes in chemokine and cytokine expression in muscle were most likely part of the adaptive mechanisms contributing to the stress response (Fig. 3). IL-16 is a lymphocyte chemoattractant factor also classified as an “alarmin” (Rider et al., 2017). IL-6 and IL-8 are also regarded as myokines released from muscle in response to contractions (Brandt & Pedersen, 2010). Muscle-derived IL-6 may mediate some of the anti-inflammatory and insulin-sensitizing effects of physical exercise (Covarrubias Anthony & Horng, 2014).

Another striking result of our study is the upregulation of transcripts related to the carbohydrate metabolic pathway, e.g., transcripts encoding PDK4 (an inactivator of pyruvate dehydrogenase complex; targeted by 1 cis-transcriptional module in the LT and 3 modules in the ST), PFKFB3 (a glycolysis regulator; targeted by 1 cis-transcriptional module in the LT and 2 cis-transcriptional modules in the ST) and SLC25A25 (a mitochondrial ATP transporter; targeted by 2 cis-transcriptional modules in the LT). This illustrates a switch in energy metabolism in the muscles of animals exposed to exercise and psychological stress towards anaerobic metabolism to support ATP production for muscle contraction. PDK4 plays a pivotal role in controlling metabolic flexibility (Zhang et al., 2014), and its expression is increased in response to moderate intensity exercise.

Analysis of molecular networks also highlighted the contribution of the response to oxygen levels/hypoxia in the response to stress, albeit via different transcripts and different contractile and metabolic muscle types. Consistently, the TF hypoxia inducible factor (HIF1A) was identified as a main common regulator of the DEGs. This could be a signature of oxygen imbalance or of the physical activity imposed on the cows. Thus, it may not be surprising that the expression of PFKB3, a downstream target of HIF, was upregulated. Hypoxia was also demonstrated to cross-talk with the Notch signalling pathways, which regulate satellite cell quiescence and self-renewal (Liu et al., 2012). Since quiescent satellite cells have a low metabolic rate, fewer mitochondria and anaerobic metabolism, this is likely part of the adaptive signature of muscle to stress. Thus, the combined signatures of hypoxia, the Notch signalling pathway (Fukada et al., 2007), MYOD1 downregulation (Kopan, Nye & Weintraub, 1994), and IRFD1 (an inducer of regenerative myogenesis) downregulation further indicate that quiescent satellite cells are stress targets and most likely physical activity targets.

Putative effects on meat quality

Finally, the transcriptomic muscle response to preslaughter stress may have an impact on meat quality through energy metabolism and hypoxia. Indeed, anaerobic glycolysis is highly relevant to beef quality since it is involved in postmortem protein degradation and hence beef tenderization during meat ageing (reviewed by Maltin et al., 2003). This process is regulated by the decline in muscle pH due to the conversion of glycogen into lactate following the lack of oxygen after slaughtering. Stress was shown to markedly affect meat tenderness by increasing postmortem ultimate pH (Purchas, 1990) due to the depletion of glycogen stores by stress prior to slaughtering, which leads to dark-cutting meats. Reliable indicators of the occurrence of high pH and preslaughter stress were identified in the sarcoplasmic proteome of muscle (Fuente-Garcia et al., 2019; Sentandreu et al., 2021). They were mainly involved in metabolism, chaperone- and stress-related processes, muscle contractility/fibre organization, and transport activities. In our study, several genes encoded by the DEGs and the common regulators or targets of the DEGs were located in bovine QTLs of the meat and carcass group associated with meat quality parameters known to be impacted by stress: muscle pH (DLL4, ERBB2), shear force (CDIPT, CEBPD, DNAJB4, GPAM, RAB3IL1, MAPK1, and TLR4), tenderness score (ADRB2 CDIPT RAB3IL1, and IFNG), compression (ADRB2), and juiciness (ATP1B1, DFFB, and RAB3IL1). Some genes were also mapped in QTLs associated with marbling score and with lipid class contents. However, the relationships between transcript levels and these meat quality parameters remain to be studied.