Transcriptomic analysis of skeletal muscle from beef cattle exposed to illicit schedules containing dexamethasone: identification of new candidate biomarkers and their validation using samples from a field monitoring trial

Growth promoters (GPs) such as the glucocorticoid dexamethasone (DEX) and the β-adrenergic agonist clenbuterol (CLEN) are still used abusively in beef cattle production. Transcriptomic markers for indirect detection of such GPs have been discussed in either experimentally treated animals or commercial samples separately. In the present study we examine the transcriptomic signature of DEX alone or in combination with CLEN in skeletal muscle of experimentally treated beef cattle, and, furthermore, compare them with previously screened commercial samples from a field-monitoring study, as well as with proteomics data representing the same set of samples. Using DNA microarray technology, transcriptomic profiling was performed on 12 samples representing three groups of animals: DEX (0.75 mg/animal/day, n = 4), a combination of DEX (0.66 mg/animal/day) and CLEN (from 2 to 6 mg/animal/day, n = 4) and a control group (n = 4). Analyses showed the differential expression of 198 and 39 transcripts in DEX and DEX–CLEN groups, respectively. Both groups had no common modulated genes in between, neither with the proteomics data. Sixteen candidate genes were validated via qPCR. They showed high correlation with the corresponding microarray data. Principal component analysis (PCA) on both the qPCR and normalised microarray data resulted in the separation of treated animals from the untreated ones. Interestingly, all the PCA plots grouped the DEX-positive samples (experimental or commercial) apart from each other. In brief, this study provides some interesting glucocorticoid-responsive biomarkers whose expression was contradictory to what is reported in human studies. Additionally, this study points out the transcriptomic signature dissimilarity between commercial and experimentally treated animals.


Introduction
The misuse of growth promoters (GPs) such as anabolic hormones, corticosteroids and β 2 -agonists in cattle is a major concern for food safety because of the public health interest. The synthetic glucocorticoid dexamethasone (DEX) is known to be illicitly used in meat cattle either alone or in combination with other active principles, such as the β 2 -agonist clenbuterol (CLEN). It is known that low doses of DEX interfere with endogenous cortisol synthesis and metabolism, therefore it results in improved feed intake, increased weight gain, reduced feed conversion ratio, reduced nitrogen retention, and increased water retention and fat content (Möstl et al. 1999;Courtheyn et al. 2002;König et al. 2006;Cannizzo et al. 2013). As regards β 2 -agonists, they elicit their growth-promoting effect by enhancing protein synthesis and cell hypertrophy, by inhibition of proteolysis in the muscle tissue and induction of lipolysis in the adipose tissue, which is known as the 'repartitioning effect'. These effects may result in a reduction of carcass fat by up to 40% and an increase of carcass protein content by up to 40%, yielding a consistent advantage for the meat industry (Courtheyn et al. 2002;Leporati et al. 2014). Pharmacologically, the repartitioning and induced muscle growth of CLEN are mediated via interaction with β-adrenergic receptors (β-AR; Rothwell & Stock 1987); however, this effect soon becomes attenuated due to decreased density (Huang et al. 2000) or desensitisation of the β-AR in skeletal muscle upon long-term exposure to β 2 -agonists (Badino et al. 2005). On the other hand, glucocorticoids are thought to reverse the homologous down-regulation of the β-AR number and mRNA expression (Mak et al. 1995), and hence are considered as part of the strategy to enhance the anabolic effects of β-adrenergic agonists (Huang et al. 2000). The effect of continuous administration of low doses of DEX and CLEN on animals' muscle mass and performance has already been illustrated (Odore et al. 2006;Biancotto et al. 2013).
The use of various illicit schedules, such as newly designed drugs or cocktails containing lower GP concentrations, makes the mission of control authorities more difficult, especially when these low GP concentrations bypass the threshold limits of current official detection methods (Cantiello et al. 2007;. For that reason, looking for alternative detection methods was inevitable. In recent years, some pilot studies investigated the effect of DEX, alone or in combination with CLEN, on different biological parameters, in order to highlight direct or indirect markers that could be used for screening purposes. In this context, regulation of tissue βadrenergic and glucocorticoid receptors in veal calves has been studied following repeated exposure to DEX or CLEN (Odore et al. 2007). Moreover, the effect of a DEX/CLEN-containing cocktail on serum immunoglobulins, lymphocyte proliferation and cytokine gene expression in veal calves has also been evaluated (Cantiello et al. 2007). Others have used protein expression changes (Stella et al. 2011), thymus morphology and serum cortisol concentration (Vascellari et al. 2012), or corticosteroid profiling of urine using coupled LC-MS/MS (Biancotto et al. 2013), as indirect markers to detect illegally administered GPs. On the other hand, commercial beef samples were also used to look for diagnostic signatures for DEX treatment in a field monitoring study (Pegolo et al. 2012). Collectively, high-throughput-omic methodologies, such as epigenomics, genomics, transcriptomics (e.g. the whole-transcriptome expression profile using DNA microarrays), proteomics and metabolomics, have recently been incorporated in this field of research, similar to other biological sciences (Riedmaier & Pfaffl 2013).
The purpose of this study was threefold: (1) to measure the changes in cattle skeletal muscle transcriptome induced by treatment with DEX alone, or in combination with CLEN, similar to those illegally performed in the field; (2) to identify a set of differentially expressed genes (DEGs) that could be used as an indirect biological marker for illicit treatment abuse; and (3) to understand, through the comparison with transcriptomic data from field monitoring samples (Pegolo et al. 2012) and proteomic results (Stella et al. 2011), if the single-approach strategy (i.e., proteomics, transcriptomics or analytical chemistry) could be enough to define a universal panel of biomarkers, or if a multi-approach strategy is needed. Finally, another question of this study was to discover whether experimentally treated animals respond differently from those present in the field.

Animals and experimental design
Twenty-four clinically healthy Charolais bulls (18-20 months old) were used. Animals were weighed, housed in ventilated stables and all the experimental procedures were carried out according to European Union animal welfare legislation. The experiment began after 3 weeks of acclimatisation. The animals were randomly divided into three groups of eight animals each. The first group was used as a control (CTR); the second group was treated with DEX, administered via feed 0.75 mg per capita for 42 days (group DEX); the third one (DEX-CLEN) was administered DEX via feed (0.66 mg per capita for 21 days) in combination with an increasing dose of CLEN, e.g. 2 mg per capita during the first week, 4 mg per capita in the second week, and 6 mg per capita during the third and fourth weeks (28 days in total; Figure 1). The products to be administered were dissolved in water, and the desired dosage was achieved by mixing 15 ml of water containing an appropriate concentration of each drug with 100 g of feed. This feed was carefully offered to each animal, ensuring that no significant residue remained in the feeder. In addition, feed conversion index (FCI) and body weight were recorded all over the study.

Sample collection and RNA extraction
At the slaughterhouse, small biceps brachii muscle tissue specimens were sampled from all animals. Muscle samples were immediately frozen in vessels containing liquid nitrogen (within 1 min of removal) and stored at -80°C prior to subsequent analyses. Total RNA was isolated by TRIzol ® reagent (Life Technologies, Grand Island, NY, USA) and subsequently purified using the RNeasy Mini kit (Qiagen, Milan, Italy), according to the manufacturer's instructions. To avoid genomic DNA contaminations, on-column DNase digestion with the RNase-free DNase set (Qiagen) was performed. Total RNA concentration was determined using the NanoDrop ND-1000 UV-Vis spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and its quality was measured by the 2100 Bioanalyzer and RNA 6000 Nano kit (Agilent Technologies, Palo Alto, CA, USA). High-quality input (hybridised RNA) is essential to have an unbiased and reproducible output in terms of gene expression data; therefore, the isolated RNAs were tested for proper concentration and integrity. The best four samples of each group (total number = 12) were selected for the microarray analyses, using a specific RNA quality parameter as a criterion, e.g. RNA concentration ≥ 40 ng µl -1 Figure 1. Experimental design (treatment). Charolais beef cattle were treated via feed with dexamethasone (DEX) alone or in combination with clenbuterol (CLEN). Animals were administered DEX at a dose rate of 0.75 mg per capita for 6 weeks. One group (DEX-CLEN, n = 8), besides DEX, was administered with CLEN increasing dosages, e.g. 2 mg per capita during the first week, 4 mg per capita during the second week, and 6 mg per capita during the third and the fourth weeks (4 weeks in total). A third group served as a control (CTR). and RNA integrity number (RIN) ≥ 6.5. The mean RIN value of these 12 samples was 6.93 ± 0.51, indicating intact RNA. Additionally, tissue specimens of four illicitly DEXtreated animals (were classified as positive by LC-MS and thymus histological analyses in the study performed by Pegolo et al. 2012) were provided and fresh RNA was isolated from them, then tested for quantity and quality as before. On the bases of having positively DEX treated samples from the field, those four samples were analysed by qPCR, then implemented in further statistical analyses (see statistics and principal component analyses) along with the experimental samples.

RNA amplification, labelling and hybridisation
Sample amplification, labelling and hybridisation were performed following the Agilent One-Color Microarray-Based Gene Expression Analysis protocol. Briefly, for each individual sample 50 ng of total RNA were linearly amplified and labelled with Cy3-dCTP using Agilent Low Input Quick Amp Labeling kit. A mixture of 10 different viral polyadenylated RNAs (Spike-In Mix, Agilent) was added to each RNA sample before amplification and labelling. A purification step was applied to the labelled cRNA using the RNeasy Mini Kit (Qiagen), and sample concentration and specific activity (pmol Cy3 µg -1 cRNA) were measured. A total of 1.65 μg of labelled cRNA was fragmented by using the Gene Expression Hybridisation kit (Agilent) according to the manufacturer's instructions, and finally diluted by the addition of 55 µl of 2X GE Hybridisation buffer. A volume of 100 µl of hybridisation solution was then dispensed in the gasket slide and assembled to the microarray slide, with each slide containing four arrays. Bovine-specific oligo-arrays (Bovine V1, 4x44k G2519F, Design ID 015354; Agilent) were used. Slides were first incubated for 17 h at 65°C in a Hybridisation Oven (Agilent), then washed using wash buffers 1 and 2 according to the manufacturer's instructions. Hybridised slides were scanned at 5 µm resolution using a G2565BA DNA microarray scanner (Agilent). Default settings were modified in order to scan the same slide twice at two different sensitivity levels (XDR Hi 100% and XDR Lo 10%). A general workflow of the microarray experiment is reported in Figure 2. Microarray data have been deposited in the NCBI's Gene Expression Omnibus (GEO) and are accessible through the GEO Series accession number GSE61934 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE61934).

Normalisation of microarray data
Data were extracted and the background was subtracted using the default settings of the Agilent's Feature Extraction Software version 9.5.1. Extracted data were normalised and processed as previously described by Giantin et al. (2014). A further filtering step was carried out by removing probes that reported missing values or no reactivity (flag equal to 0) in at least 50% of samples. Missing values (probes with Feature Extraction flag equal to 0) were imputed by using the microarray data analysis tool TIGR Multiple Array Viewer (TMEV; Saeed et al. 2003). In addition, raw microarray data (30 samples) from a monitoring study were downloaded from GEO (GSE26318; Pegolo et al. 2012). Experimental and monitoring samples (n = 42) were grouped, normalised and filtered together in one run to avoid any possible analytical bias, and then used in subsequent analyses.

Data mining and pathway analyses
The functional analysis of DEGs list was performed using the Ingenuity Pathway Analysis (IPA) online platform (see http://www.ingenuity.com; Qiagen).

Quantitative real time PCR
Sixteen target genes (Table 1) and three housekeeping (HK) genes (RPLP0, GAPDH and TBP) were chosen for the external validation of microarray findings by qPCR. To increase the number of data points in the statistical and correlation analyses between DNA microarray and qPCR, samples from all 24 animals in the study were included in the qPCR analysis. Gene-specific primers (see Table S1 in the Supplemental data online), encompassing one intron, and the most appropriate Universal Probe Library (UPL) probe were designed by using the UPL Assay Design Centre web service (Roche Applied Science, Indianapolis, IN, USA). Firststrand cDNA was synthesised from 0.15 µg of total RNA using the High Capacity cDNA Reverse Transcription Kit (Life Technologies) according to the manufacturer's protocol and stored at -20°C until further use. Overall, qPCR reactions (10 µl final volume) consisted of 1X LightCycler 480 Probe Master (Roche Applied Science), 300 or 600 nM forward and reverse primers (Integrated DNA Technology, Italy) according to the assay set-up, 200 nM human UPL probe (final concentrations) and 2.5 µl of 1:7.5 diluted cDNA (15 ng µl -1 ). Each qPCR analysis was performed, in duplicate, in a LightCycler 480 Instrument (Roche Applied Science) using the standard PCR conditions (an activation step at 95°C for 10 min; 45 cycles at 95°C for 10 s and at 60°C for 30 s; a cooling step at 40°C for 30 s) and LightCycler 480 clear plates (Roche Applied Science). To determine the efficiency of each qPCR assay, non-template and no reverse transcription controls were included on each plate. Moreover, standard curves obtained by amplifying eight threefold serial dilution of the same cDNA pool were used. Data were analysed with the LightCycler480 software release 1.5.0 (Roche Applied Science) using the second derivative or fit-point method. mRNA relative quantification was performed by the ΔΔCq method (Livak & Schmittgen 2001).

Statistics and principal component analysis (PCA)
To identify DEGs, a two-class unpaired test was implemented in the program SAM (Significance Analysis of Microarrays) release 4.0 (Tusher et al. 2001), enforcing a false discovery rate (FDR) of 5% with a fold change (FC) threshold of 2. All other statistical tests (linear regression, non-parametric Spearman correlation analysis and Mann-Whitney test) were carried out by the GraphPad Prism 5 software (San Diego, CA, USA). Statistical significance was set at p < 0.05.
Using TMEV, a PCA analysis was carried out on our microarray processed samples (n = 12) in order to test the efficiency of the identified DEGs in differentiating treated from untreated samples. Hence, only normalised and  To visualise the multivariate response of the selected classifiers to the treatment and to see if the chosen set of gene markers could result in the best separation of treated and untreated animals, another dynamic PCA approach, based on qPCR results (relative quantification values, RQs), was performed. The objective was to test if the proposed set of the 16 genes could differentiate DEX-treated samples (experimental or commercial) from the untreated ones. For that purpose, the later PCA was performed using RQ values corresponding to the CTR (n = 8) and DEX (n = 8) groups, along with the four samples claimed to be positive for DEX from the monitoring study. The PCA on qPCR RQ values was executed by GenEx v.5 software (Bergkvist et al. 2010), adopting the following settings: mean centre scaling, Ward's algorithm and Manhattan distance. Finally, another PCA analysis was carried out on our microarray processed (n = 12) and all the monitoring samples (n = 30) or the DEX-positive ones (n = 4) in order to test if they could group into different clusters based on their DNA microarray raw data. For each PCA, samples were grouped together in one set, and a blind PCA was carried out.

Animal health status and growth performance
The health status of all the experimental animals was satisfactory throughout the experiment. Except for DEX and/or CLEN, no other drugs were administered during the experimental procedure. Information about the animals' performance and feed conversion index have been presented in detail by Biancotto et al. (2013).

Microarray quality control and data analyses
The comparison of the DEX group with the CTR one (DEX versus CTR) resulted in a list of 198 downregulated transcripts, representing 123 characterised transcripts and 75 estimated sequence tags (ESTs: see Table S2 in the Supplemental data online). Estimated sequence tag transcripts were excluded from further analyses due to limited or unavailable information on annotation. Among the 123 down-regulated genes, 29 had FC > 4, and the highest FC (-14.77-fold) was noticed for the LOC509808 (LIPG) gene.
Following the comparison between DEX-CLEN and CTR groups (DEX-CLEN versus CTR), the analysis with SAM resulted in a much shorter list of DEGs. A total of 39 DEGs, representing 21 characterised transcripts and 18 ESTs (thereby excluded from further analyses), was identified. Among the 21 genes, 16 and five were respectively up-and down-regulated (see Table S3 in the Supplemental data online). There were no overlapping DEGs between the two comparisons.
To describe better the transcriptome functional modifications and to mine through the obtained DEG lists, IPA was used to find out and explore pathways, networks and biofunctions somewhat modulated by the used GPs. Data analysis identified some networks and biofunctions, e.g. molecular and cellular functions, physiological system development and canonical pathways (see Table S4 in the Supplemental data online for a detailed IPA report). Moreover, glucocorticoid receptor and xenobiotic metabolism signalling pathways were highlighted among the canonical pathways' output following the IPA core analysis performed on our DEGs (Figure 3).
To reduce the complexity of a long list of significant genes, an additional approach was to focus on those transcripts that seemed more relevant to our study. Through an extensive literature screening, using the NCBI database and gene summaries present in GeneCards of the human gene database (see http://www.genecards.org), together with some predictions suggested by IPA, we collected information about the direct function of each gene and the corresponding target tissue (when available). This strategy helped to refine the DEG lists and to group genes into three main categories, e.g. glucocorticoids responsive, skeletal muscle related and coagulation cascade-related genes (see Table S5 in the Supplemental data online). A set of 16 genes shown to be modulated in DEX (15 genes: CYP1A1, FKBP5, CCL24, RASD1, MYOC, PFKFB4, MEDAG, SULT1A1, LIPG, CRISPLD2, C1QA, FGL2, CCDC80, C7 and MMP2) and DEX-CLEN groups (one gene: HSPA8) were considered for the final validation analysis by qPCR.

Confirmatory qPCR analysis
To cross-validate our platform performance and identify some potential biomarkers, a relative quantification approach by using qPCR was a must. All the aforementioned 16 genes were found to be significantly regulated in the muscle tissue, comparing treated with control animals (Figure 4). A Spearman rank correlation test was then performed for each target gene, comparing qPCR RQ values with the corresponding microarray's intensities ( Figure 5). In the DEX group, 14 out of 15 genes showed high correlation coefficients (Spearman rho > 0.8; p < 0.001), while only one gene (MEDAG) was significantly correlated, but with a lower correlation coefficient (rho = 0.685; p = 0.01). In the DEX-CLEN group, a significant relationship was only found for HSPA8 (rho = 0.734; p < 0.01; Table 1). The correlation made on the average FC obtained with DNA microarray and qPCR technologies showed a high and significant correlation coefficient (Spearman's r = 0.754, p < 0.0001), with a slope of linear regression of 0.94 ( Figure 6).

Clustering and PCA
The first two components, which accounted for 81.124% of the total variance, clearly identified and distinguished DEX from CTR samples, while DEX-CLEN samples were distributed along the xand y-axes with no clear grouping trend ( Figure 7A).
By using the RQ values of the proposed 16 gene markers, the PCA was able to distinguish three groups of samples, e.g. CTR (n = 8), DEX-treated animals (n = 8) and the four DEX-positive samples from the monitoring study. The first two principal components of greatest variation covered 94.60% of the total variance. Quite unexpectedly, samples from DEX group and DEX-positive samples from the monitoring (MO) study did not group together on the PCA plot ( Figure 7B). Likewise, samples clustered independently from each other also on the relative dendrogram ( Figure 8).
To understand better why DEX-treated animals clustered apart from each other, and particularly to verify whether this behaviour was a treatment-specific response or a technical bias, a broad dynamic PCA, using the raw microarray data of all the samples (present study: n = 12, field monitoring: n = 30), was performed ( Figure 9A). In this case, the first two components, in compliance with Pegolo et al. (2012), accounted for 66.5% of the total variance and separated the monitoring samples (MO-1 and MO-2) from our experimental samples (EXP). The y-axis, representing 24.3% of the total variance, separated our 12 samples from group MO-2, including the four DEX-positive samples and other unknown samples.
To confirm these findings further and to avoid any possible distortion in the PCA plot due to many unknown samples in the monitoring study, the PCA was repeated Figure 3. (colour online) Ingenuity pathways analysis (IPA). Canonical pathways tree following a core analysis by IPA and using human, mouse and rat databases as a background. The glucocorticoid receptor signalling pathway is highlighted. between our 12 samples (EXP) and the four DEX-positive samples (MO-POS), using the differentially regulated gene list (11 484 unique transcripts) following a oneclass SAM analysis, and enforcing an FDR of 0%. Again, we found that our samples and the monitoring ones were grouped apart from each other ( Figure 9B).

Discussion
Although the development of additional detection methods for GPs' abuse in beef cattle has already been faced by different 'omic' disciplines, such as genomics, transcriptomics, proteomics and metabolomics, the 'gold standard' detection technology or the perfect biomarkers panel is seemingly out of reach for the moment. In the present study we investigated the effect of some GPs on the muscle's transcriptome of beef cattle. The foremost objective was to detect possible changes in gene expression and to compare these transcriptional effects with those obtained applying a proteomic approach (Stella et al. 2011) and MS-based analytical investigations (Biancotto et al. 2013). Moreover, another goal was to test our suggested transcriptional biomarkers against prevalidated DEX-positive samples (Pegolo et al. 2012). Shortly, a set of potential transcriptional biomarkers was identified and cross-validated with an independent method. However, the comparison of our DEGs list with that of the proteomics study revealed no common gene-coded protein or pathway(s) in between. Finally, after plotting our samples on PCA against DEX- Figure 4. Validation of 16 differentially expressed target genes. Expression levels (relative quantification (RQ) values) of genes, relative to non-treated (CTR) samples determined by qPCR. Significance was tested by using the one-way analysis of variance (ANOVA) (Kruskal-Wallis test) followed by Dunn's post-hoc test. Multiple bars show means ± SE. *, **, ***Statistically significant differences of DEX versus CTR values (p < 0.05, 0.01, 0.001, respectively). positive samples coming from a monitoring study, our proposed target genes did not group the experimental and monitoring samples together; on the contrary, they were clustered separately. These distinct results are hereby discussed more in depth.
Indirect biological markers have a comparable cost, higher output and high sensitivity compared with other methods (Balizs & Hewitt 2003;Carraro et al. 2009). Based on results here obtained, the use of illicit schedules containing DEX alone could be reliably identified, with high confidence, by using 15 genes; on the other hand, the use of a DEX-CLEN combination was eventually identified by the HSPA8 gene only. The identification of 198 DEGs in the DEX-treated group and only 39 ones in the DEX-CLEN group could suggest a CLEN-masking effect upon DEX. Counteraction between DEX and CLEN, where the latter caused mild attenuation of the effects of DEX on some physiological parameters, has already been reported (Huang et al. 2000). Several genes were broadly down-regulated in the DEX-treated group, and we will shortly discuss only the ones showing the highest response in terms of fold changes and/or relevance to the purpose of the study. The most down-regulated gene (-15.32-fold) is the endothelial lipase (LIPG) gene. Endothelial lipase is a member of the triglycerides (TG) lipase gene family, showing a significant phospholipase activity on high-density lipoprotein (HDL) particles (Goldberg 1996;Jaye et al. 1999;Griffon et al. 2006); furthermore, LIPG gene inactivation has been shown to increase HDL levels (Qiu et al. 2007). Although the relationship between LIPG and lipid metabolism is well documented, there is not much information in the literature regarding any relationship between this gene and glucocorticoids. However DEX has been shown to decrease LIPG activity in the liver of rats (Peinado-Onsurbe et al. 1991). Moreover, a recent study reported that DEX administration in horses decreased the expression of genes involved in hormone signalling, cholesterol synthesis and steroidogenesis (Ing et al. 2014). The well-studied and evident effect of DEX on lipid metabolism (Zhou & Cidlowski 2005;Xu et al. 2009;Martin et al. 2009;Campbell et al. 2011) could be a reason to hypothesise a relationship between LIPG and DEX treatment.
The cytochrome P450 1A1 gene was down-regulated by 10.39-fold. This gene is involved in oxidative drug metabolism (Beresford 1993) and it depends on the aryl hydrocarbon receptor (AhR) for its regulation. It has been reported that CYP1A1 expression in adult human hepatocytes was negatively regulated by DEX at the protein level, but no effects were noticed upon mRNA (Monostory et al. 2005). On the other hand, DEX was shown to have no inhibitory potency on the CYP1A1 level either in human hepatocytes (Vrzal et al. 2009) or in rainbow trout (Burkina et al. 2013). In the same context, DEX was recently proved to suppress CYP1A1 transactivation in gene reporter assays (Stejskalova et al. 2013).
FKBP5 (also known as FKBP51) can act as an important determinant of the responses to steroids, especially to glucocorticoids in stress and mood disorders in humans (Kang et al. 2008;Jääskeläinen et al. 2011). Up-regulation of FKBP5 in response to corticosteroid use has been consistently demonstrated in many studies Almon et al. 2005;Tissing et al. 2007), and has been associated with the loss of efficacy of corticosteroids (Fisher et al. 2005). Surprisingly, in the present study FKBP5 was down-regulated (-7.14-fold). Most of the information about FKBP5 is closely related to humans, and very little information is available about this gene in animals, if any. However, it has been recently reported that FKBP5 is down-regulated following 21-day treatment with the corticosteroid prednisolone in a collagen-induced arthritis mouse model (Ellero-Simatos et al. 2014). It should be emphasised that illicit schedules in cattle substantially differ, in terms of dosage and duration of administration, from those used in humans, and that the glucocorticoid signalling system is highly stochastic and differs greatly between tissues (Kino 2007).
Dexamethasone-induced Ras-related protein 1 (RASD1) is a member of the Ras family of proteins that is usually activated following the administration of corticosteroids (Kemppainen & Behrend 1998;Tu & Wu 1999). Surprisingly, we observed a down-regulation of this gene (-3.05-fold). To our knowledge, there is no available information about RASD1 gene expression in cattle. In humans, RASD1 mRNA is constitutively expressed in many tissues such as brain, heart, liver and kidney (Kemppainen & Behrend 1998;Tu & Wu 1999;Fang et al. 2000), while no records are available for skeletal muscle. The inactivation of RASD1 and its correlation with resistance to DEX, as a consequence of methylation, was recently discussed by Nojima et al. (2009). A thorough explanation of these contradictory findings needs further investigations and comparative studies to be performed in order to understand if there are species-and/or tissue-related variations.
The chemokine (C-C motif) ligand 24 (CCL24) is predominantly involved in chemotaxis of T-cells compatible with an anti-inflammatory effect of glucocorticoids (Ehrchen et al. 2007). In few studies performed on cattle, CCL24 was shown to be normally down-regulated in the gestation period of female cows (Oliveira et al. 2010;Laporta et al. 2014). However, data from human and experimental animal models, in accordance with our findings, showed that CCL24 was down-regulated following DEX treatment (Goleva et al. 2008;Luesink & Jansen 2010;Louten et al. 2012). Overall, these findings are in concordance with the down-regulation (-4.84-fold) we observed following DEX treatment.
The myocilin, trabecular meshwork inducible glucocorticoid response (MYOC) gene was shown to be moderately down-regulated (-2.36-fold). This result is apparently in contrast with other studies where MYOC was up-regulated upon DEX treatment in both human and cattle anterior segments of the eye, namely the trabecular meshwork cells (Taniguchi et al. 2000;Ishibashi et al. 2002;Rozsa et al. 2006;Figure 6. Overall comparison between microarray and qPCR data, expressed as estimated fold-changes, referring to the final 16 candidate transcript biomarkers. ***p < 0.001. Mao et al. 2011). Nevertheless, it was recently reported (Kumar et al. 2013) that MYOC showed decreased expression in trabecular meshwork cells isolated from eyes of mice treated with steroids. Furthermore, the effect of DEX upon MYOC gene expression varies greatly from one tissue matrix to another (Morgan et al. 2014). One of the other interesting results was the downregulation (-3.44-fold) of the cysteine-rich secretory protein LCCL domain containing 2 (CRISPLD2) gene. This transcript has only recently been considered as a glucocorticoid-responsive gene (Himes et al. 2014). What is interesting here is that CRISPLD2 showed an over-expression pattern in some primary human cell lines following 24-48 h of DEX treatment (Masuno et al. 2011;Greer et al. 2013;Himes et al. 2014), while it was down-regulated in our study. Very little information is actually available for this gene in the literature, except for humans. The present contradictory result might be due to the different illicit protocol (GP combination) and/or the DEX dosages here used.
Likewise to CRISPLD2, the complement component 7 (C7), a key component of the adaptive immunity (Actor et al. 2001;Mashruwala et al. 2011), has been recently considered as a glucocorticoid-responsive gene (Himes et al. 2014). In our experimental conditions, C7 was found to be down-regulated (-2.12-fold). This could be explained by the fact that glucocorticoids are known to repress the expression of adaptive immune-related genes Franchimont 2004;Azuma 2010). The only up-regulated gene in our validated set of genes was the heat shock proteins A8 (HSPA8), which has been examined as a glucocorticoid-induced gene in other model systems (Smoak & Cidlowski 2006). However, no information relating this gene to CLEN is present in the literature so far.
A second purpose of the present work was the comparison between our DEG lists (transcriptomics data) with the one referring to differentially expressed proteins (e.g., proteomics data), for which the same muscle samples from the same animal were used (Stella et al. 2011). Interestingly, the comparison resulted into no common differentially regulated genes/protein in between. In cattle, the proteomic approach has been explored as a tool for the detection of protein expression patterns in skeletal muscles (Keady et al. 2013;Stella et al. 2014), body fluids (Draisci et al. 2007;Della Donna et al. 2009;Guglielmetti et al. 2014) and some target organs such as liver (Gardini et al. 2006). Until recently, there was an implicit assumption in the systems biology literature about the existence of a proportional relationship between mRNA and protein expressions measured from tissue. However, analysis of mRNA and protein expression data from the same cells under similar conditions have failed to show a high correlation between the two domains in multiple studies (Pascal et al. 2008;Ghazalpour et al. 2011). Moreover, small non-coding RNAs (miRNA) and post-translational modifications such as phosphorylation, SUMOylation and ubiquitination have been shown to modulate the expression, regulation, stability and function of glucocorticoid molecular targets and pathways (Duma et al. 2006;De Iudicibus et al. Figure 8. Dendrogram showing the hierarchical clustering of the three experimental groups; controls (CTR, B4:1-B4:8), animals administered with dexamethasone (DEX, B4:9-B4:16) and monitoring DEX-positive cattle (MO: S21, S54, S56 and S57). This tree represents the similarity between genes and/or samples based on the gene expression profiles (RQ values) measured by qPCR assays. Average linkage and Euclidean distance were used as a clustering method and distance measure, respectively. 2013). Therefore, finding no single gene overlap was somewhat surprising, but further confirms that the transcriptomic and proteomic approaches are independent of each other, even if the target tissue is the same. This agrees with Timperio et al. (2009), who reported that proteomics and transcriptomics data seldom overlapped when compared upon analysis of liver samples from two different Bos taurus breeds. Hence, the need for an integrated approach against GPs' misuse in cattle seems to be required to improve the effectiveness of the indirect biomarker approach for screening purposes.
On the microarray level, the defined DEGs were able to distinguish our experimental groups one from another. This has an effect on the efficiency of the analyses and the robustness of the obtained data. In addition, the validated and proposed set of biomarkers (16 genes, whose level of expression was measured by qPCR) grouped the field monitoring DEX-positive samples away from our experimental DEX group on the PCA; furthermore, they were also clustered into different batches on a dendrogram. This unexpected behaviour of the monitoring samples prompted us to decide to use the whole microarray raw dataset of all the monitoring samples along with ours to check whether or not monitored commercial animals still behave differently from the experimentally treated ones. Interestingly, all our samples grouped together and distinctively away from the other monitoring samples, which were represented on the PCA as Figure 9. PCA of the bovine skeletal muscle gene expression profiles. (A) PCA plot showing the grouping of the 12 experimental samples distinctively away from commercial field monitoring samples (MO-1 and MO-2). The plot was created following a one-class SAM analysis (via TMEV) on normalised microarray intensities by using a differentially regulated gene list (11 484 unique transcripts) and a false discovery rate (FDR) of 0%. The first two components accounted for 42.2% and 24.3%, respectively. (B) PCA plot excluding all the monitoring study samples except for the four dexamethasone (DEX)-positive ones; the xand y-axes cover 60.8% and 8.8% of the total variance, respectively. The ellipse on the right side represent the 12 experimental samples, while the four samples on the left side of the plot are the four DEX-positive ones from the monitoring study. previously reported by Pegolo et al. (2012). Moreover, this result was confirmed after repeating the PCA using our 12 samples against only the four DEX-positive samples. Despite the efforts made to decrease any possible outliers affecting samples' distribution on the PCA, we again had the two groups apart. Once again, variables such as breed (perhaps), diet and breeding conditions (more probably), the use of different DEX-containing illicit protocols and/or new cocktails containing different GPs could be the reason(s) for this conflict. Indeed, intraspecies transcriptomic comparison is needed to reveal the differences in drug response between animals under field and experimental conditions.

Conclusions
The present work has shown that despite the attractions of comparative profiling of transcripts and proteins on a global 'omic' scale, there are obvious biological and technical differences preventing transcriptomics and proteomics from having a convergence. Moreover, the idea of having a 'universal set of biomarkers' that can be applied to all illicitly treated cattle still seems elusive at the moment. Indeed, one 'stand-alone' technology does not suffice to gain a comprehensive understanding of the biological system complexity. An approach that incorporates the various 'omic' platforms and their data would be the key to solve this puzzle.