BACKGROUND
in the european union the use of growth promoting substances in livestock production is prohibited according ec directive 96/ <dig>  <cit> .
to ensure compliance with this legislation, requirements for monitoring are described in ec directive 96/ <dig>  <cit> .
at national level, legislations are implemented in residue monitoring programs regulating sampling of animal matrices and residue analysis therein to guarantee fair trade, food safety and public health.
residue analysis in livestock production is in general based on chemical  <cit> , immunochemical or biological  <cit>  screening methods followed by mass spectrometry based confirmation methods.
although this strategy seems to work for synthetic anabolic steroids, problems arise when compounds that also occur naturally are used.
abuse of naturally occurring hormones is hard to prove since most of these substances are strongly metabolized in vivo.
moreover, metabolites are not always known or are present in levels not significantly different from highly fluctuating endogenous levels.
this makes it difficult to prove fraudulent use based on quantification of natural occurring compounds.
nowadays, it is observed that misuse of growth promoters in cattle fattening moves towards these natural steroids and steroid esters.
moreover, inspections of livestock farms in the netherlands occasionally result in the finding of feed or herbal additives and preparations containing so-called prohormones.
prohormones are compounds that exhibit limited or no hormonal action by themselves, however they are direct precursors of active hormones and indirectly affect natural hormone levels.
dehydroepiandrosterone  is such a prohormone and is the most abundant occurring precursor of both androgens and estrogens in humans  <cit> .
it is claimed that orally taken dhea improves muscle strength and is therefore illicitly used in sports to enhance performance and appearance  <cit> .
looking for alternatives to support evidence of illegal use of growth promoting substances, gene expression analysis can be an attractive new approach.
several studies demonstrated changes in mrna expression in bovine tissues upon treatment with growth promoters after performing real-time rt-pcr analysis on a limited number of preselected genes  <cit> .
untargeted transcriptomics approaches using microarrays allow gene expression analysis of thousands of genes simultaneously as well as identification of  biomarkers for screening  <cit> .
moreover, microarray data can provide mechanistic insights in cellular processes and pathways and can be used for classification of compounds with the same mode of action   <cit> .
comparative microarray analysis is therefore in potential a promising screening tool for growth promoter abuse and in particular for prohormones of which the mode of action in cattle is sometimes unclear.
in recent work we used a metabolomics approach to compare urine profiles of control and dhea exposed bovines  <cit> .
this revealed several urinary steroid phase i and phase ii metabolites which are potential biomarkers for dhea treatment.
in the present study we investigated the feasibility of monitoring prohormone abuse at the mrna level using liver tissue from the same animal experiment.
gene expression profiles of control and dhea treated animals were compared to determine differentially expressed genes and to identify biomarkers for dhea treatment.
methods
animals and treatment
male frisian bull calves were purchased at the local market and housed for 2- <dig> weeks before the start of the experiment.
treatment with dhea was repeated in three independent experiments using identical treatment and sampling schedules.
each of the three experiments consisted of two animals of which one was orally  treated with capsules containing  <dig> mg dhea  and the other was injected intramuscularly  with  <dig> mg dhea dissolved in  <dig> ml miglyol  <dig> .
administrations were performed seven times, at 24-hour intervals.
im treated animals  were 9- <dig> months old and po treated animals  8- <dig> months old.
control animals were included in the first  and third  experiment.
an overview of the experimental setup, age and weights of the bovines is shown in table  <dig>  twenty-four hours after the last treatment, the animals were sacrificed and liver tissue was collected, snap-frozen in liquid nitrogen and stored at -80°c until use.
the experimental work was approved by the animal ethics committee of ghent university, belgium, in accordance with local ethical requirements.
microarray analysis
total rna was extracted from tissues by homogenization in trizol  and mixed with chloroform.
the lysate was centrifuged at  <dig> × g for  <dig> minutes at 4°c and the aqueous phase was transferred to be mixed with isopropanol which precipitates total rna.
after centrifuging  the pellet was washed with 75% ethanol and resuspended in rnase free water.
upon extraction the rna was purified according to the rneasy mini kit protocol .
after purification, rna integrity was determined spectroscopically  and by gel electrophoresis.
only rna with a260/ <dig> and a260/ <dig> ratios above  <dig>  was used for amplification.
to generate fluorescently-labeled crna, the agilent low rna input fluorescent linear amplification kit  was used according to the manufacturer's protocol.
in short,  <dig> μg of total rna was reverse transcribed using t <dig> tagged oligo-dt primer and labeled with cy <dig> or cy <dig> .
liver rnas of the treated and control animals were individually labeled with cy <dig> and rna of all  <dig> control animals was pooled and labeled with cy <dig>  after purification with the rneasy mini kit , label efficiency and yield were determined using a nanodrop spectrophotometer .
a mixture of  <dig> μg of cy3-labeled crna and  <dig> μg of cy5-labeled crna was hybridized onto a 44k bovine oligo microarray , using agilent's gene expression hybridization kit.
hybridization was performed at 65°c for  <dig> hours in a hybridization oven with rotation function .
upon hybridization, microarrays were washed and dried according to the agilent's instructions and fluorescence measurements were performed using an agilent technologies g2565b microarray scanner.
data processing
fluorescence intensities were quantified using feature extraction  <dig>  software .
data were imported in genemaths xt  <dig>   and signals below two times background were excluded from further analysis.
subsequently, the data was normalized as described by pellis et al.
<cit> .
this normalization included correction for the random error, with the median cy <dig> signal for each individual spot.
secondly, correction for the systematic error was performed with the median value of the overall cy <dig> signal.
after normalization, principal component analysis  was performed to visualize differences between groups and t-test statistics were performed to test for differential expression.
microarray data was floored by adjusting low intensity spots to a threshold value of  <dig>  hereby reducing the number of less reliable genes.
next, spot intensities were 2log transformed and each gene was mean centered versus all samples.
based on these 2log transformed data differentially regulated genes were selected with a p-value < <dig>  and a fold change > <dig>   in each of the three treatment replicates versus the average from the control animals.
hierarchical clustering of the differentially regulated genes was performed using cluster and treeview software  <cit> .
raw microarray data of the present study have been submitted to arrayexpress  and are stored under experiment accession number a-mexp- <dig>
gene set enrichment analysis
gene set enrichment analysis  is a tool to identify and analyse the differential expression of biologically relevant sets of genes that share common biological functions  <cit> .
using gsea, the differentially regulated genes observed in dhea treated animals versus one control group  were validated by evaluating this gene set by comparing the same dhea treated animals versus the other control group .
therefore, separate gene sets were generated of the differentially expressed genes of respectively im treated animals  and os treated animals  versus the controls of experiment  <dig> as well as the controls of experiment  <dig>  for example, the transcripts found significantly up-regulated when comparing dhea im versus control group  <dig> were included in the gene set "dhea_im_vs_ctr1_up".
in a similar way other gene sets were created for up- as well as down-regulated genes.
next, gsea ranks all genes on the microarray on differential expression between dhea exposed and controls using signal to noise statistics, resulting in a list with up-regulated genes at the top and down-regulated genes at the lower end of the list.
each of the pre-defined gene sets was tested against this list and gsea calculated whether the genes in the gene set are randomly distributed, enriched at the top or at the lower end of the ranked list.
permutations were performed on gene sets and gene sets were considered significantly affected when the p-value was below  <dig>  and the false discovery rate  below  <dig> , accordingly to gsea recommendations  <cit> .
RESULTS
principal component analysis  and selection of differentially regulated genes
in the present study the potential strengths as well as the pitfalls of microarray experiments using calves from real-life practice were investigated.
three small animal experiments were performed independently using bull calves purchased at the local market.
in this way the experimental setup was taking into account the inherent variability needed to investigate the usefulness of bovine-specific microarrays as a screening tool for prohormone abuse in veterinary control.
for obvious ethical reasons larger numbers of bovines treated with banned substances could not be justified.
upon microarray hybridization and data normalization, unsupervised principal component analysis  was performed to visualize differences between liver profiles of control and treated animals.
figure  <dig> shows the pca-plot which is based on the three largest components, representing  <dig> % of the total variance.
although there is variation in gene expression profiles of livers of animals treated with dhea, they are clearly discriminated from those of the controls.
however, large differences are observed between the two control groups, whereas the exposed animals  and the control animals of the first experiment  are separated along the x-axis while the control bovines of experiment three  and the exposed bovines are mainly separated along the z-axis.
based upon the outcome of this pca, further analysis was focused on comparison of the im and po treated animals versus either the total control population as well as the two control groups separately.
differentially regulated genes were selected using t-test statistics.
a p-value < <dig>  and a difference of at least  <dig>   fold change, versus the control average, observed in all three biological treatment replicates  were used as criteria for the selection of differentially expressed genes.
an overview of the differentially regulated genes found in the im and po treated animals is shown in figure  <dig>  a total of  <dig> and  <dig> genes were found to be regulated in im and po treated animals as compared to the total control group, respectively.
only one of these genes  was found differentially expressed  in im as well as po treated animals.
a hierarchical cluster diagram of all differentially regulated genes is presented in figure  <dig>  since many probes were spotted twice or more on the microarray the  <dig> and  <dig> genes found regulated are represented by  <dig> and  <dig> spots respectively.
a detailed description of all regulated genes is listed in additional file  <dig>
of the  <dig> differentially expressed genes in response to the im dhea treatment,  <dig> represented unidentified transcripts and  <dig> were encoding for proteins with poorly known or unknown function .
among the  <dig> transcripts that encode for known proteins,  <dig> transcripts are involved in immune response and inflammatory processes.
of these latter transcripts xbp <dig>  mx <dig>  lbp, serpina <dig>  ccl <dig>  card <dig> and pigr were found up-regulated and ankrd <dig>  lyz and dmbt were down-regulated.
the remaining transcripts are involved in various processes like cell growth and proliferation , formation of tight-junctions , tumor suppression , cell proliferation and cell adhesion  intra-cellular signalling  and cell-cell interactions .
regarding metabolism, the gstp <dig> gene was found to be down regulated > <dig>  fold in all im treated animals.
gstp <dig> mediates glutathione conjugation and plays an important role in detoxification of xenobiotics as well as in uptake and transport of numerous hydrophobic endogenous compounds like steroids  <cit>  moreover, it has been observed in mouse that the gstp <dig> gene contains androgen receptor binding sites which regulate gstp <dig> activity in response to androgens  <cit> .
comparison of po treated animals versus all control animals revealed a total of  <dig> differentially expressed genes of which  <dig> represent unidentified transcripts or encoded for proteins with an unknown function.
again a substantial number of the differentially regulated genes are involved in immune response of which lilra <dig>  thbs, clec6a and fut <dig> were found up-regulated and ccl <dig> and dbmt were down-regulated.
other differentially regulated genes are involved in peptidase inhibition , g-protein signalling  and amino acid transport .
also regulated is the short/branched chain acyl-coa dehydrogenase  gene, a member of the acyl-coa dehydrogenase enzyme family which is involved in fatty acid metabolism.
this could point towards regulation of fatty acid metabolism and is supported by a study in which dhea administration to rats showed significant regulation of genes involved in fatty acid metabolism, including the very long chain acyl-coa gene which is also a member of the acyl-coa dehydrogenase enzyme family  <cit> .
overall it can be stated that the majority of regulated genes are involved in immune response for both po as well as im treated animals which is in line with numerous studies reporting the significant immune modulatory properties of dhea  <cit> .
in principle the above listed genes are potential biomarkers for dhea treatment.
on the one hand, we are aware of the small number of animals used in this study which hampers proper statistics and substantially increases the chance of missing dhea-responsive genes or detecting false-positive genes.
on the other hand, combining and comparing the data of three independently performed experiments will limit the risk of false-positive genes considerably and results in identification of only the most robustly regulated genes.
therefore, we assessed whether the genes differentially expressed in animals treated with dhea via one administration route versus animals of one control group would also be affected when compared 1) to other control animals and 2) by the other administration route.
to deal with these issues we applied the statistics of gene set enrichment analysis .
in this way statistical power could be improved and regulated gene sets were tested for their robustness.
gene set enrichment analysis
for gsea we used the genes found to be differentially expressed when the exposed animals are compared with the two control groups separately.
for dhea im,  <dig> and  <dig> genes were differentially regulated versus ctr <dig> and ctr <dig>  respectively, whereas for dhea po,  <dig> and  <dig> genes were regulated.
as shown in figure  <dig>  a relatively small number of genes was affected per treatment in both comparisons versus ctr <dig> and comparisons versus ctr <dig>  apparently, only a low proportion of genes showed a significant up-regulation or down-regulation of  <dig>  or more in both comparisons.
gsea was used as a tool to discriminate dhea treated animals from non-treated animals on the basis of gene sets generated from genes found to be differentially regulated .
gene sets were compared to the whole experimental dataset and gsea calculated whether genes within a gene set are randomly distributed, enriched at the top or at the bottom of the ranked list  <cit> .
the advantage of this gsea approach is that no cut-off is used for determination of differentially regulated genes.
using the whole experimental data set makes that alterations are viewed for as a group of genes instead for individual genes.
gene sets can be significantly affected while changes in expression of individual genes are relatively subtle.
for example, the transcripts found to be significantly up-regulated when comparing dhea im versus control group  <dig>  were included in the gene set 'dhea_im_vs_ctr3_up'.
gsea analysis, using this gene set, showed that the genes where highly enriched in dhea im treated animals when comparing versus the ctr1-group .
as shown in figure 4c, most of the genes are distinctly up-regulated in dhea im treated animals, although also variation in gene expression of the individual animals is observed.
in a similar way, the other gene sets were compared versus the other 'independent' control group and results are summarized in table  <dig>  this cross validation showed that  <dig> out of  <dig> gene sets were significantly enriched  when dhea treated animals were compared versus an 'independent' control group.
moreover, gene sets generated on the basis of dhea im treated animals showed significant enrichment in dhea po treated animals and vice versa.
in total  <dig> out of  <dig> gene set comparisons showed significant enrichment, suggesting an overlap in gene expression profiles from im and po treated animals which most likely include genes that are differentially expressed irrespective of the manner of dhea administration.
* significantly regulated gene set.
although in vivo transcriptomics data of dhea in liver is limited  <cit> , depreter et al.
identified  <dig> genes which were found to be up-regulated in rat liver  <cit> .
gsea analysis showed significant enrichment of this gene set in dhea im and po treated animals when these were compared with the controls of experiment  <dig> .
these results illustrate that gsea is a powerful approach for comparative analysis of gene expression data obtained in different settings.
controlled experiments with bovines have resulted in the identification of genomics based biomarkers which potentially can be used for screening for hormones  <cit> .
however, when examining bovines from real-life practice, one is dealing with biological variation like age, genetic background, environment, nutrition and disease history.
in the current study, this biological variation was deliberately included and was mainly reflected by the large differences in gene expression profiles of the control populations tested.
the two control groups in this study showed substantial age differences i.e.
the animals in the ctr <dig>  ctr <dig> and dhea-treated groups are  <dig> months, 13- <dig> months and 8- <dig> months in age, respectively.
nevertheless, for the dhea im and po treated animals, sets of respectively  <dig> and  <dig> genes were found differentially expressed when compared to all controls using standard statistics.
these two groups of genes are specific for im and po treatment, respectively, and independent from biological factors like age.
however, gsea results showed a correlation between gene expression profiles of im and po treated animals, suggesting that there are also effects independent from the route of administration.
this is in line with our earlier performed metabolomics study showing large similarities in urine metabolite profiles of im and po treated animals as well as metabolites specific for the route of administration  <cit> .
hence, for application of transcriptomics based screening of bovines for hormones in practice, the treatment effect should be filtered out from differences in gene expression due to inherent biological variation.
here it was shown that microarray gene expression profiling in combination with statistical methods like gsea are able to distinguish gene expression profiles of dhea-treated animals from non-treated control animals.
it should be noted that this experiment comprised small numbers of animals and follow up experiments are required to gain statistical power and to obtain a better description of dhea specific gene sets.
furthermore, the behaviour of such a growth promoter specific gene set should be studied in a broad spectrum of untreated control animals from daily practice, to assure the robustness of the gene set.
this underlines the need to obtain more liver gene expression profiles of control animals from slaughterhouses.
CONCLUSIONS
the present study showed that identification of genomic biomarkers for dhea treatment in cattle is hampered by the large biological variability as compared to genomics experiments with inbred strains of rodents under well-defined laboratory conditions.
however, gene expression profiling using whole genome microarrays in combination with predefined gene sets and statistical methods like gsea showed to be a promising approach to screen for hormone abuse in livestock production.
for application in practice however, a better genomic description of the control population as well as growth promoter specific gene set are needed.
authors' contributions
jcwr was involved in the setting up the experiments, carrying out the microarray experiments and writing of the draft manuscript.
aacmp and mjg participated in setting up the animal study and in writing the draft manuscript.
pjmh participated in data analysis and contributed to the draft manuscript.
jmvh carried out the animal experiments and helped with sampling.
mwfn initiated the study, participated in its design and helped to draft the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
differentially expressed genes.
tables of differentially expressed genes
click here for file
acknowledgements
this project was financially supported by the dutch ministry of agriculture, nature and food quality.

 in the european union the use of growth promoting substances in livestock production is prohibited according ec directive 96/. to ensure compliance with this legislation, requirements for monitoring are described in ec directive 96/. at national level, legislations are implemented in residue monitoring programs regulating sampling of animal matrices and residue analysis therein to guarantee fair trade, food safety and public health. residue analysis in livestock production is in general based on chemical, immunochemical or biological screening methods followed by mass spectrometry based confirmation methods. although this strategy seems to work for synthetic anabolic steroids, problems arise when compounds that also occur naturally are used. abuse of naturally occurring hormones is hard to prove since most of these substances are strongly metabolized in vivo. moreover, metabolites are not always known or are present in levels not significantly different from highly fluctuating endogenous levels. this makes it difficult to prove fraudulent use based on quantification of natural occurring compounds. nowadays, it is observed that misuse of growth promoters in cattle fattening moves towards these natural steroids and steroid esters. moreover, inspections of livestock farms in the netherlands occasionally result in the finding of feed or herbal additives and preparations containing so-called prohormones. prohormones are compounds that exhibit limited or no hormonal action by themselves, however they are direct precursors of active hormones and indirectly affect natural hormone levels. dehydroepiandrosterone is such a prohormone and is the most abundant occurring precursor of both androgens and estrogens in humans. it is claimed that orally taken dhea improves muscle strength and is therefore illicitly used in sports to enhance performance and appearance. looking for alternatives to support evidence of illegal use of growth promoting substances, gene expression analysis can be an attractive new approach. several studies demonstrated changes in mrna expression in bovine tissues upon treatment with growth promoters after performing real-time rt-pcr analysis on a limited number of preselected genes. untargeted transcriptomics approaches using microarrays allow gene expression analysis of thousands of genes simultaneously as well as identification of biomarkers for screening. moreover, microarray data can provide mechanistic insights in cellular processes and pathways and can be used for classification of compounds with the same mode of action. comparative microarray analysis is therefore in potential a promising screening tool for growth promoter abuse and in particular for prohormones of which the mode of action in cattle is sometimes unclear. in recent work we used a metabolomics approach to compare urine profiles of control and dhea exposed bovines. this revealed several urinary steroid phase i and phase ii metabolites which are potential biomarkers for dhea treatment. in the present study we investigated the feasibility of monitoring prohormone abuse at the mrna level using liver tissue from the same animal experiment. gene expression profiles of control and dhea treated animals were compared to determine differentially expressed genes and to identify biomarkers for dhea treatment. methods animals and treatment male frisian bull calves were purchased at the local market and housed for 2- weeks before the start of the experiment. treatment with dhea was repeated in three independent experiments using identical treatment and sampling schedules. each of the three experiments consisted of two animals of which one was orally treated with capsules containing mg dhea and the other was injected intramuscularly with mg dhea dissolved in ml miglyol. administrations were performed seven times, at 24-hour intervals. im treated animals were 9- months old and po treated animals 8- months old. control animals were included in the first and third experiment. an overview of the experimental setup, age and weights of the bovines is shown in table twenty-four hours after the last treatment, the animals were sacrificed and liver tissue was collected, snap-frozen in liquid nitrogen and stored at -80°c until use. the experimental work was approved by the animal ethics committee of ghent university, belgium, in accordance with local ethical requirements. microarray analysis total rna was extracted from tissues by homogenization in trizol and mixed with chloroform. the lysate was centrifuged at × g for minutes at 4°c and the aqueous phase was transferred to be mixed with isopropanol which precipitates total rna. after centrifuging the pellet was washed with 75% ethanol and resuspended in rnase free water. upon extraction the rna was purified according to the rneasy mini kit protocol. after purification, rna integrity was determined spectroscopically and by gel electrophoresis. only rna with a260/ and a260/ ratios above was used for amplification. to generate fluorescently-labeled crna, the agilent low rna input fluorescent linear amplification kit was used according to the manufacturer's protocol. in short, μg of total rna was reverse transcribed using t tagged oligo-dt primer and labeled with cy or cy. liver rnas of the treated and control animals were individually labeled with cy and rna of all control animals was pooled and labeled with cy after purification with the rneasy mini kit, label efficiency and yield were determined using a nanodrop spectrophotometer. a mixture of μg of cy3-labeled crna and μg of cy5-labeled crna was hybridized onto a 44k bovine oligo microarray, using agilent's gene expression hybridization kit. hybridization was performed at 65°c for hours in a hybridization oven with rotation function. upon hybridization, microarrays were washed and dried according to the agilent's instructions and fluorescence measurements were performed using an agilent technologies g2565b microarray scanner. data processing fluorescence intensities were quantified using feature extraction software. data were imported in genemaths xt and signals below two times background were excluded from further analysis. subsequently, the data was normalized as described by pellis et al.. this normalization included correction for the random error, with the median cy signal for each individual spot. secondly, correction for the systematic error was performed with the median value of the overall cy signal. after normalization, principal component analysis was performed to visualize differences between groups and t-test statistics were performed to test for differential expression. microarray data was floored by adjusting low intensity spots to a threshold value of hereby reducing the number of less reliable genes. next, spot intensities were 2log transformed and each gene was mean centered versus all samples. based on these 2log transformed data differentially regulated genes were selected with a p-value < and a fold change > in each of the three treatment replicates versus the average from the control animals. hierarchical clustering of the differentially regulated genes was performed using cluster and treeview software. raw microarray data of the present study have been submitted to arrayexpress and are stored under experiment accession number a-mexp- gene set enrichment analysis gene set enrichment analysis is a tool to identify and analyse the differential expression of biologically relevant sets of genes that share common biological functions. using gsea, the differentially regulated genes observed in dhea treated animals versus one control group were validated by evaluating this gene set by comparing the same dhea treated animals versus the other control group. therefore, separate gene sets were generated of the differentially expressed genes of respectively im treated animals and os treated animals versus the controls of experiment as well as the controls of experiment for example, the transcripts found significantly up-regulated when comparing dhea im versus control group were included in the gene set "dhea_im_vs_ctr1_up". in a similar way other gene sets were created for up- as well as down-regulated genes. next, gsea ranks all genes on the microarray on differential expression between dhea exposed and controls using signal to noise statistics, resulting in a list with up-regulated genes at the top and down-regulated genes at the lower end of the list. each of the pre-defined gene sets was tested against this list and gsea calculated whether the genes in the gene set are randomly distributed, enriched at the top or at the lower end of the ranked list. permutations were performed on gene sets and gene sets were considered significantly affected when the p-value was below and the false discovery rate below, accordingly to gsea recommendations.  principal component analysis and selection of differentially regulated genes in the present study the potential strengths as well as the pitfalls of microarray experiments using calves from real-life practice were investigated. three small animal experiments were performed independently using bull calves purchased at the local market. in this way the experimental setup was taking into account the inherent variability needed to investigate the usefulness of bovine-specific microarrays as a screening tool for prohormone abuse in veterinary control. for obvious ethical reasons larger numbers of bovines treated with banned substances could not be justified. upon microarray hybridization and data normalization, unsupervised principal component analysis was performed to visualize differences between liver profiles of control and treated animals. figure shows the pca-plot which is based on the three largest components, representing % of the total variance. although there is variation in gene expression profiles of livers of animals treated with dhea, they are clearly discriminated from those of the controls. however, large differences are observed between the two control groups, whereas the exposed animals and the control animals of the first experiment are separated along the x-axis while the control bovines of experiment three and the exposed bovines are mainly separated along the z-axis. based upon the outcome of this pca, further analysis was focused on comparison of the im and po treated animals versus either the total control population as well as the two control groups separately. differentially regulated genes were selected using t-test statistics. a p-value < and a difference of at least fold change, versus the control average, observed in all three biological treatment replicates were used as criteria for the selection of differentially expressed genes. an overview of the differentially regulated genes found in the im and po treated animals is shown in figure a total of and genes were found to be regulated in im and po treated animals as compared to the total control group, respectively. only one of these genes was found differentially expressed in im as well as po treated animals. a hierarchical cluster diagram of all differentially regulated genes is presented in figure since many probes were spotted twice or more on the microarray the and genes found regulated are represented by and spots respectively. a detailed description of all regulated genes is listed in additional file of the differentially expressed genes in response to the im dhea treatment, represented unidentified transcripts and were encoding for proteins with poorly known or unknown function. among the transcripts that encode for known proteins, transcripts are involved in immune response and inflammatory processes. of these latter transcripts xbp mx lbp, serpina ccl card and pigr were found up-regulated and ankrd lyz and dmbt were down-regulated. the remaining transcripts are involved in various processes like cell growth and proliferation, formation of tight-junctions, tumor suppression, cell proliferation and cell adhesion intra-cellular signalling and cell-cell interactions. regarding metabolism, the gstp gene was found to be down regulated > fold in all im treated animals. gstp mediates glutathione conjugation and plays an important role in detoxification of xenobiotics as well as in uptake and transport of numerous hydrophobic endogenous compounds like steroids moreover, it has been observed in mouse that the gstp gene contains androgen receptor binding sites which regulate gstp activity in response to androgens. comparison of po treated animals versus all control animals revealed a total of differentially expressed genes of which represent unidentified transcripts or encoded for proteins with an unknown function. again a substantial number of the differentially regulated genes are involved in immune response of which lilra thbs, clec6a and fut were found up-regulated and ccl and dbmt were down-regulated. other differentially regulated genes are involved in peptidase inhibition, g-protein signalling and amino acid transport. also regulated is the short/branched chain acyl-coa dehydrogenase gene, a member of the acyl-coa dehydrogenase enzyme family which is involved in fatty acid metabolism. this could point towards regulation of fatty acid metabolism and is supported by a study in which dhea administration to rats showed significant regulation of genes involved in fatty acid metabolism, including the very long chain acyl-coa gene which is also a member of the acyl-coa dehydrogenase enzyme family. overall it can be stated that the majority of regulated genes are involved in immune response for both po as well as im treated animals which is in line with numerous studies reporting the significant immune modulatory properties of dhea. in principle the above listed genes are potential biomarkers for dhea treatment. on the one hand, we are aware of the small number of animals used in this study which hampers proper statistics and substantially increases the chance of missing dhea-responsive genes or detecting false-positive genes. on the other hand, combining and comparing the data of three independently performed experiments will limit the risk of false-positive genes considerably and results in identification of only the most robustly regulated genes. therefore, we assessed whether the genes differentially expressed in animals treated with dhea via one administration route versus animals of one control group would also be affected when compared 1) to other control animals and 2) by the other administration route. to deal with these issues we applied the statistics of gene set enrichment analysis. in this way statistical power could be improved and regulated gene sets were tested for their robustness. gene set enrichment analysis for gsea we used the genes found to be differentially expressed when the exposed animals are compared with the two control groups separately. for dhea im, and genes were differentially regulated versus ctr and ctr respectively, whereas for dhea po, and genes were regulated. as shown in figure a relatively small number of genes was affected per treatment in both comparisons versus ctr and comparisons versus ctr apparently, only a low proportion of genes showed a significant up-regulation or down-regulation of or more in both comparisons. gsea was used as a tool to discriminate dhea treated animals from non-treated animals on the basis of gene sets generated from genes found to be differentially regulated. gene sets were compared to the whole experimental dataset and gsea calculated whether genes within a gene set are randomly distributed, enriched at the top or at the bottom of the ranked list. the advantage of this gsea approach is that no cut-off is used for determination of differentially regulated genes. using the whole experimental data set makes that alterations are viewed for as a group of genes instead for individual genes. gene sets can be significantly affected while changes in expression of individual genes are relatively subtle. for example, the transcripts found to be significantly up-regulated when comparing dhea im versus control group were included in the gene set'dhea_im_vs_ctr3_up'. gsea analysis, using this gene set, showed that the genes where highly enriched in dhea im treated animals when comparing versus the ctr1-group. as shown in figure 4c, most of the genes are distinctly up-regulated in dhea im treated animals, although also variation in gene expression of the individual animals is observed. in a similar way, the other gene sets were compared versus the other'independent' control group and results are summarized in table this cross validation showed that out of gene sets were significantly enriched when dhea treated animals were compared versus an'independent' control group. moreover, gene sets generated on the basis of dhea im treated animals showed significant enrichment in dhea po treated animals and vice versa. in total out of gene set comparisons showed significant enrichment, suggesting an overlap in gene expression profiles from im and po treated animals which most likely include genes that are differentially expressed irrespective of the manner of dhea administration. * significantly regulated gene set. although in vivo transcriptomics data of dhea in liver is limited, depreter et al. identified genes which were found to be up-regulated in rat liver. gsea analysis showed significant enrichment of this gene set in dhea im and po treated animals when these were compared with the controls of experiment. these results illustrate that gsea is a powerful approach for comparative analysis of gene expression data obtained in different settings. controlled experiments with bovines have resulted in the identification of genomics based biomarkers which potentially can be used for screening for hormones. however, when examining bovines from real-life practice, one is dealing with biological variation like age, genetic background, environment, nutrition and disease history. in the current study, this biological variation was deliberately included and was mainly reflected by the large differences in gene expression profiles of the control populations tested. the two control groups in this study showed substantial age differences i.e. the animals in the ctr ctr and dhea-treated groups are months, 13- months and 8- months in age, respectively. nevertheless, for the dhea im and po treated animals, sets of respectively and genes were found differentially expressed when compared to all controls using standard statistics. these two groups of genes are specific for im and po treatment, respectively, and independent from biological factors like age. however, gsea results showed a correlation between gene expression profiles of im and po treated animals, suggesting that there are also effects independent from the route of administration. this is in line with our earlier performed metabolomics study showing large similarities in urine metabolite profiles of im and po treated animals as well as metabolites specific for the route of administration. hence, for application of transcriptomics based screening of bovines for hormones in practice, the treatment effect should be filtered out from differences in gene expression due to inherent biological variation. here it was shown that microarray gene expression profiling in combination with statistical methods like gsea are able to distinguish gene expression profiles of dhea-treated animals from non-treated control animals. it should be noted that this experiment comprised small numbers of animals and follow up experiments are required to gain statistical power and to obtain a better description of dhea specific gene sets. furthermore, the behaviour of such a growth promoter specific gene set should be studied in a broad spectrum of untreated control animals from daily practice, to assure the robustness of the gene set. this underlines the need to obtain more liver gene expression profiles of control animals from slaughterhouses.  the present study showed that identification of genomic biomarkers for dhea treatment in cattle is hampered by the large biological variability as compared to genomics experiments with inbred strains of rodents under well-defined laboratory conditions. however, gene expression profiling using whole genome microarrays in combination with predefined gene sets and statistical methods like gsea showed to be a promising approach to screen for hormone abuse in livestock production. for application in practice however, a better genomic description of the control population as well as growth promoter specific gene set are needed. authors' contributions jcwr was involved in the setting up the experiments, carrying out the microarray experiments and writing of the draft manuscript. aacmp and mjg participated in setting up the animal study and in writing the draft manuscript. pjmh participated in data analysis and contributed to the draft manuscript. jmvh carried out the animal experiments and helped with sampling. mwfn initiated the study, participated in its design and helped to draft the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 differentially expressed genes. tables of differentially expressed genes click here for file acknowledgements this project was financially supported by the dutch ministry of agriculture, nature and food quality.
changes in gene expression profiles were studied in livers of bull calves treated orally or intramuscularly with mg dhea versus two control groups, using bovine 44k dna microarrays. as determined by principal component analysis, large differences in liver gene expression profiles were observed between treated and control animals as well as between the two control groups. when comparing the gene expression profiles of po and im treated animals to that of all control animals, the number of significantly regulated genes was and respectively. however, a better insight in the genomic variability of the control population is a prerequisite in order to define growth promoter specific gene sets that can be used as robust biomarkers in daily practice. this study showed that identification and application of genomic biomarkers for screening of hormone abuse in livestock production is substantially hampered by biological variation. on the other hand, it is demonstrated that comparison of pre-defined gene sets versus the whole genome expression profile of an animal allows to distinguish dhea treatment effects from variations in gene expression due to inherent biological variation. therefore, dna-microarray expression profiling together with statistical tools like gsea represent a promising approach to screen for hormone abuse in livestock production. in contrast to controlled genomics studies, this work involved bovines purchased at the local market on three different occasions with ages ranging from to months, thereby reflecting the real life inter-animal variability due to differences in age, individual physiology, season and diet. this cross validation, showed that out of the gene sets were significantly enriched in dhea treated animals when compared to an'independent' control group. illegal use of natural prohormones like dehydroepiandrosterone is hard to prove since prohormones are strongly metabolized in vivo. within the european union the use of growth promoting agents in animal production is prohibited. in the present study, we investigated the feasibility of a novel effect-based approach for monitoring abuse of dhea. for im and po treated calves, gene sets were generated of genes that were significantly regulated compared to one control group and validated versus the other control group using gene set enrichment analysis.
 within the european union the use of growth promoting agents in animal production is prohibited. illegal use of natural prohormones like dehydroepiandrosterone is hard to prove since prohormones are strongly metabolized in vivo. in the present study, we investigated the feasibility of a novel effect-based approach for monitoring abuse of dhea. changes in gene expression profiles were studied in livers of bull calves treated orally or intramuscularly with mg dhea versus two control groups, using bovine 44k dna microarrays. in contrast to controlled genomics studies, this work involved bovines purchased at the local market on three different occasions with ages ranging from to months, thereby reflecting the real life inter-animal variability due to differences in age, individual physiology, season and diet.  as determined by principal component analysis, large differences in liver gene expression profiles were observed between treated and control animals as well as between the two control groups. when comparing the gene expression profiles of po and im treated animals to that of all control animals, the number of significantly regulated genes was and respectively. for im and po treated calves, gene sets were generated of genes that were significantly regulated compared to one control group and validated versus the other control group using gene set enrichment analysis. this cross validation, showed that out of the gene sets were significantly enriched in dhea treated animals when compared to an'independent' control group.  this study showed that identification and application of genomic biomarkers for screening of hormone abuse in livestock production is substantially hampered by biological variation. on the other hand, it is demonstrated that comparison of pre-defined gene sets versus the whole genome expression profile of an animal allows to distinguish dhea treatment effects from variations in gene expression due to inherent biological variation. therefore, dna-microarray expression profiling together with statistical tools like gsea represent a promising approach to screen for hormone abuse in livestock production. however, a better insight in the genomic variability of the control population is a prerequisite in order to define growth promoter specific gene sets that can be used as robust biomarkers in daily practice. 