BACKGROUND
farm animals have been highly selected for favorable production traits such as lean growth rate, feed efficiency, and prolificacy in pigs.
in parallel their robustness has tended to decrease, as shown by the sensitivity to diseases, locomotor weakness or behavioral problems  <cit> .
the hypothalamic-pituitary-adrenocortical  axis plays a major role in physiological homeostasis including metabolism, brain function and behavior, the immune system and inflammatory processes.
together with the autonomic nervous system, it is also involved in stress and adaptive responses to environmental challenges.
on the basis of available literature, we have hypothesized that the hpa axis could contribute to the trade-off between production and robustness traits, and that genetic variation in hpa axis activity could be used to select more robust animals  <cit> .
hpa axis activity shows a large variation among individuals and genetic influences are well documented  <cit> .
for example, in pigs, the sensitivity of the adrenal cortex to adrenocorticotropic hormone  and the production of corticosteroid binding globulin , the carrier of cortisol in blood are the two main mechanisms responsible for genetic differences in circulating cortisol levels  <cit> .
in a previous paper, hazard et al.
<cit>  have studied at the gene expression level the molecular mechanisms of genetic differences in the adrenal gland response to acth.
however little is known about the individual variation in the biological activity of cortisol, the main glucocorticoid hormone, and the genetic mechanisms involved.
corticosteroid hormones exert their biological actions via two intracellular receptors  that, upon activation by their ligand, influence the expression of a large number  of genes in a wide range of cell types  <cit> .
in pigs, glucocorticoid receptor polymorphisms have been described with their functional consequences .
the present experiment was designed to explore in pigs individual variation in the biological response to cortisol in order to identify possible biomarkers of this sensitivity.
to that end, juvenile pigs were submitted to an acth challenge to release cortisol and the kinetics of the response was measured with biological variables and with gene expression analysis in blood cells.
taken together with biological information, this approach will serve as an important step to understand hpa axis regulation and will identify key genes involved in signaling pathways relevant to stress responses.
the final goal of this work is to develop a strategy for further genetic studies in order to overcome the unfavourable consequences of stress in farm animals.
animals and methods
animals, treatment and blood sampling
all animal use was performed under european union and french legislation .
the protocol and procedures were approved by the local  ethics committee .
experimental animals were  <dig> piglets  randomly selected from  <dig> litters  of purebred large white pigs and produced in  <dig> successive batches raised  <dig> weeks apart.
they were weaned at  <dig> weeks and animals from 2– <dig> litters were mixed at weaning in each post-weaning pen.
experimental animals were not isolated from their littermates.
at  <dig> weeks, each animal was injected in the neck muscles with  <dig> μg of synthetic 1– <dig> acth  diluted in  <dig> ml of  <dig>  % saline.
injections occurred from 10:00–11: <dig> am to avoid nycthemeral variations.
blood samples were collected before the injection  and  <dig> ,  <dig>  and  <dig>  hours later.
at each time, animals were slightly restrained on their back in such a way that the effect on their stress level can be regarded as insignificant.
two blood samples were then taken by puncture of a jugular vein in vacutainer Ⓡ tubes with  <dig> g needles.
the whole handling procedure lasted less than  <dig> sec.
one  <dig> ml tube with lithium heparin was used for chemical biology.
after centrifugation , plasma aliquots were frozen at – <dig> °c until analysis.
one  <dig> ml tube with edta  was used for blood cell count and an aliquot  was mixed with the same volume of dl buffer , frozen at – <dig> °c for  <dig> h and then at – <dig> °c until analysis for gene expression.
biological analyses
cortisol was measured by direct automated immunoassay .
glucose and free fatty acid , were measured by colorimetry with an abx pentra  <dig> clinical chemistry analyzer from horiba medical .
blood cell counts were measured with a ms-9- <dig> hematology analyzer from melet schloesing laboratories , calibrated for pig blood by the manufacturer.
blood cell count variables included: white cells count, proportion of lymphocytes, monocytes and granulocytes, red cells count, percentage of hematocrit, concentration of hemoglobin, red cells width and volume, concentration of platelets and platelets width and volume.
the biological variables contained thus  <dig> variables measured on the  <dig> pigs.
in addition, birth and weaning weights were also measured for each pig.
outlying observations were visually identified and treated as missing data.
missing data were imputed using a k-nearest neighbour  imputation .
to ensure normality, cortisol, platelets and white cells count were log <dig> transformed and ffa was transformed using the square-root.
batch effects were removed by aligning the within-batch medians for all measurements.
rna extraction and whole blood gene expression analysis
a total rna isolation and purification was done according to the manufacturer’s instructions using the nucleospin rna blood kit  followed by dnase treatment.
the quality of each rna sample was checked through the bioanalyser agilent  <dig>  and low-quality rna preparations were discarded .
microarray description
a porcine microarray gpl <dig>  was used to hybridize the rna samples as previously described  <cit> .
this microarray contained  <dig>  spots.
among them,  <dig> were negative controls and  <dig> were used for aligning.
one probe was duplicated twice on each array.
thus the microarray contained  <dig>  unique porcine probes.
after quality control, quantile normalization and filtering,  <dig>  transcripts were found to be expressed in blood in our experimental conditions.
hybridization protocol
blood samples from  <dig> female pigs from only  <dig> batches were used.
each of the  <dig> arrays used contained  <dig> microarrays which corresponded to the  <dig> observations of two individuals, each from one batch.
this design secured the kinetics of the response for each individual and prevented confounding effects between batch and array.
after quality control and filtering,  <dig>  probes were kept and log <dig> transformed.
technical biases were handled by aligning the within-array medians for all genes and by a quantile normalization within animal .
no missing data were reported.
normalized data were submitted to geo/ncbi with the gse <dig> accession number.
statistical analyses
all analyses were performed with the r software, version  <dig> .
<dig>  <cit> .
multilevel approach
a multilevel approach was used to investigate the relationships between the repeated measurements while taking advantage of multivariate approaches.
the multilevel approach, as described by liquet et al.
<cit> , uses a split-up variation coming from the mixed-model framework.
let x=i= <dig> …,n,t∈{ <dig> + <dig> + <dig> +24},k= <dig> …,p be the  observation matrix  on n animals with  <dig> times of measurements .
x can be split up as:
x=x..⏟offset term+xb⏟between-animal deviation+xw⏟within-animal deviation
the matrix x.. represents the offset term defined as 1nx..t where 1n is a vector of length n containing only ones and x..t= .
xb is the between-animal matrix of size  defined by concatenating 14xbit for each animal into xb with xbit=xi.1−x.. <dig> …,xi.p−x..pxi.k=14∑t∈{ <dig> + <dig> + <dig> +24}xitk.
xw=x−xi.
is the within-animal deviation matrix of size  with xi.
the matrix defined by concatenating the matrices 14xi.t for every animal into xi., with xi.t=xi.
<dig> …,xi.p.
by splitting the different parts of the variation in the data while taking into account the repeated measurements on the subjects, the multilevel step allows us to study the effect of different conditions within a subject separately from the variation between subjects.
this method is especially relevant when a high between-subject variability is observed in repeated data: multivariate approaches  can then be performed on xw to highlight the most relevant correlations between variables in the dataset, independently from individual variations.
statistical analysis of plasma metabolites and hormone
first, all variables were subjected to a one-way anova with repeated measures.
p-values were adjusted using a benjamini-hochberg  approach in order to control the false discovery rate   <cit> .
variables with an adjusted p-value  were then subjected to  <dig> paired t-tests to assess the difference between t= <dig> and t=+ <dig>  between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig>  the full list of p-values was adjusted using a bh approach.
in addition, the influence of sex on the biological variables was tested using a two-way anova with repeated measures including sex as a variable.
p-values were adjusted using a bh approach.
then, variability between individuals before the acth injection was studied by performing a pca on the observations at t= <dig> for variables responding to acth and birth and weaning weights.
the overall effect of acth over time was investigated with a multilevel pca as previously described.
cortisol levels at t=+ <dig> is the most relevant measure to assess the sensitivity of the adrenals to acth.
hence, correlations between biological variables at t∈{ <dig> + <dig> + <dig> +24} and the level of cortisol at t=+ <dig> were investigated using paired t-tests.
p-values were adjusted using a bh approach.
statistical analysis of the transcriptome
all transcripts were subjected to  <dig> paired t-tests to assess the difference in expression between t= <dig> and t=+ <dig>  between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig>  the full list of p-values was adjusted using a bonferroni approach.
as the bonferroni approach exerts a more stringent control than the bh approach, it was used to obtain a narrowed list of the most significantly differentially expressed  transcripts.
transcripts with at least one adjusted p-value < <dig>  among the three tests were kept.
correlations between the expression levels of different transcripts of the same gene were investigated to highlight genes for which at least 3/ <dig> of the duplicates were significantly de and had a correlation of at least  <dig>  between them.
the most significant transcripts per annotated gene were kept and the multilevel approach was used to extract the within-subject deviation matrix for further analysis.
a hierarchical ascending classification  was performed using the ward method with an euclidean distance matrix based on the correlations between genes.
this allowed for the characterization of clusters of genes having the same evolution over time.
significance of the difference between time measurements was assessed using  <dig> paired t-tests at the average gene level between t= <dig> and t=+ <dig>  between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig>  p-values were adjusted within the clusters using a bh approach.
integration
relations between the main stress variable, the cortisol, the biological variables and the gene expressions were studied using different approaches.
as for biological data, correlations between gene expression at t∈{ <dig> + <dig> + <dig> +24} and the level of cortisol at t=+ <dig> were investigated.
more precisely, the within-subject components of the transcriptomic data and of the cortisol level was extracted using the method described in  <cit>  and section “multilevel approach”.
then, a paired t-test was performed between the within-subject expression at t∈{ <dig> + <dig> + <dig> +24} and the within-subject level of cortisol at t=+ <dig>  the full list of p-values was globally adjusted using a bh approach.
in order to assess the significativity of change in expression after removing the contribution of cell population changes, a linear mixed model was fitted for every deg
xit=β0+β <dig> t+β2it+ui+εit  where xit is the expression of the deg for the animal number i  and the time step t , it is the lymphocytes/granulocytes ratio for the same experiment and ui is the individual random effect.
both time step  and it were supposed to have fixed effects on gene expression.
significance of the time effect in this model was checked by testing β <dig> t=β <dig>  for t∈{+ <dig> + <dig> +24} and correcting p-values by time point for multiple tests with a bh approach .
the effect of changes in white cell populations was also assessed by testing β2= <dig> and a correction for multiple tests was applied using a bh approach .
a multilevel pls, i.e., a pls performed on the within-subject components of the biological and transcriptomic data, as described in  <cit>  and section “multilevel approach”, was used to investigate the overall relationships between biological and transcriptomic data.
a sparse version of the pls  as described in  <cit>  was used to select a small subset of variables to explain every axis.
sequence annotation
each cdna sequence was compared to refseq_rna mammalian database using the ncbi blastn program .
resulting hits were sorted out according to their closeness to the pig genome, their coverage and sequence identity.
the selected cdna sequences were submitted to the hugo  gene nomenclature committee, using their refseq ids .
then, hugo gene symbols or official gene symbol were used as gene names.
functional enrichment and pathway analysis
functional enrichment was performed on each list of clustered genes identified by hac and on the list of genes of which expression is significatively explained by cell population changes in blood .
functional annotation of genes was provided by the biomart database  <cit> .
to set the statistical enrichment of a particular biological function, a fisher’s exact test was performed, using the list of genes expressed on the microarray as the reference list of genes.
resulting p-values were adjusted for multiple tests using a bh approach.
a minimum of  <dig> genes per gene ontology and a fdr< <dig>  were necessary to consider a biological function to be enriched.
a pathway is an interconnected arrangement of processes, representing the functional roles of genes in the genome.
functional integration of gene expression, i.e., identification of gene networks, was performed using the ‘gene ontology’ database amigo .
the significantly up- or down-regulated genes could be assembled into networks using ‘ingenuity pathway analysis’  under licence.
this application provides computational algorithms to identify the enriched biological pathways, functions and biological mechanisms of selected genes and proposes also enriched regulators as transcription factors.
a regulatory network could be constructed with the information provided by the option ‘upstream regulator’.
this option proposes a list of regulators known to have a significant effect on some of the targeted genes in the input list.
ingenuity also provides computational algorithms to identify and to dynamically generate significant biological networks.
networks are ranked by a score that takes into account the number of focus genes and the size of the networks.
this score ) indicates the probability for genes to be associated in the same network by chance.
a higher score means a smaller probability for genes to be observed in the same network by chance.
we chose networks displaying direct relationships between genes.
path designer  was used to improve the readability of the networks.
nodes added by ingenuity were discarded when they were not necessary to connect our genes of interest and the resulting network was merged with the regulatory network.
RESULTS
plasma cortisol, metabolites and hematologyfig.
<dig> mean evolution of the biological variables overtime.
: time measurement at which the expression of the variable is significantly different  from the expression at t= <dig>  vertical bars correspond to + and – sem at each time point
f value and fdr are for the test of the global time effect on each variable.
nr non relevant since the measure is the same at all time steps
as expected  <cit> , acth induced a strong cortisol response peaking  <dig> h after injection .
plasma cortisol levels at t=+ <dig> were  <dig> -fold higher than and highly correlated with basal levels .
they were significantly lower than baseline values at t=+ <dig>  and at t=+ <dig> , as can be expected from the feedback of cortisol on its own secretion.
plasma glucose levels showed a slight increase at t=+ <dig> after acth injection, but ffa levels showed a large increase at t=+ <dig> , with a strong variation across animals, and were almost back to basal levels at t=+ <dig>  the values measured at t=+ <dig> were correlated with basal values .
cortisol induces a weak increase in circulating glucose and also potentiates the effect of other counter-regulatory hormones  <cit>  and increases ffa levels via acute lipolysis  <cit> .
the data obtained here for clinical hematology measures are in the range of published values in pigs.
these values show large variations with age and breed among other sources .
although the total number of leucocytes was only marginally influenced by acth, large changes in leucocyte subpopulations were observed with an increase of the proportion of granulocytes and a decrease of the proportion of lymphocytes and monocytes at t=+ <dig> and t=+ <dig>  these effects are consistent with the literature  <cit>  and result from the redistribution of leucocytes between blood and tissues  <cit> .
red blood cell related variables  were decreased after injection of acth and remained low at t=+ <dig>  platelets were not influenced by injection of acth.
sex did not influence any of the variables .
between-subject variability at t=0
results of the pca at t= <dig> on variables responding to acth and birth and weaning weights are shown in fig.
<dig>  as red cells count , hematocrit  and hemoglobin  are redundant variables, decision was made to keep only rc for the pca.
no outliers are identified on the projection of the individuals on the two first dimensions of the pca.
on this pca, sex was not found to be related to the main variability  in the dataset.
the first dimension shows an opposition between the proportion of granulocytes  and the proportion of lymphocytes .
the second axis shows an opposition between cortisol  and birth and  weaning weights  with a strong opposition between these variables on the whole plan formed by the first and second axis.
the other variables were not correlated with either of the first two dimensions of the pca.
fig.
<dig> pca on the biological variables identified as responding to acth, the birth and weaning weights at t= <dig>  colors symbolize the sex: black = female; red = male a projection of the individuals on dimensions  <dig> and 2; b projection of the variables on dimensions  <dig> and 2; bw = birth weight; ww = weaning weight; lympho = lymphocyte ratio; mono = monocytes ratio; granulo = granulocyte ratio; rc = red cell count; gluc = glucose; ffa = free fatty acids
overall effect of the injection of exogenous acth on clinical biology variables
extraction of the within-subject data matrix prior to the application of a pca analysis allows for the separation of the observations according to their time of measurement .
the first component of the multi-level pca opposes the observations at t= <dig> and t=+ <dig>  to the observations at t=+ <dig> , this time step corresponding to the peak of cortisol.
the second component opposes the observations at t=+ <dig>  to the observations at t=+ <dig> .
the representation of the variables shows that the first axis is mainly driven by an opposition between the proportion of granulocytes, ffa, cortisol and red cell count , on one side, and lymphocytes and monocytes , on the other side.
the second axis shows an opposition between glucose  and cortisol, ffa and red count .
fig.
<dig> multilevel pca on the biological variables responding to acth.
colors symbolize the time of measurement; black: t=0; red: t=+1; green: t=+4; blue: t=+24; a projection of the individuals on dimensions 1–2; b projection of the variables on dimensions 1–2; lympho = lymphocyte ratio; mono = monocyte ratio; granulo = granulocyte ratio; rc = red cell counts; gluc = glucose; ffa = free fatty acids
specific links to the level of cortisol at t=+1
correlations between cortisol at t=+ <dig> and other variables at t= <dig> allows for the identification of variables which baseline levels may be directly or indirectly linked to the intensity of the cortisol level in response to acth, a measure of individual variation in hpa axis activity.
correlations are shown in table  <dig>  only glucose and ffa levels at t= <dig> were significantly positively correlated with the level of cortisol at t=+ <dig> .
correlations between cortisol at t=+ <dig> and other variables at t=+ <dig> t=+ <dig> or t=+ <dig> allows for the identification of variables which are directly linked to the level of cortisol when it reaches its peak during the stress response.
ffa at t=+ <dig> were positively correlated with cortisol at t=+ <dig> and glucose at t=+ <dig> t=+ <dig> and t=+ <dig> was negatively correlated with cortisol at t=+ <dig> .
no other variable was found to be significantly linked to the intensity of the cortisol level in response to acth.
<dig>
<dig>
<dig>
– <dig>
– <dig>
– <dig>
se: standard error of the correlation coefficient; in bold: significantly ≠ <dig>
differentially expressed genes
we used a comprehensive gene expression profiling by means of microarray analysis to identify clusters of genes differentially expressed in peripheral blood cells, taking into consideration the kinetics of the response with  <dig> time points .
differential analysis revealed  <dig> de transcripts  matching  <dig> unique genes .
among them,  <dig> genes were differentially expressed at t=+ <dig> ,  <dig> were differentially expressed at t=+ <dig>  and  <dig> were differentially expressed at t=+ <dig> .
the only gene de at both t=+ <dig> and t=+ <dig> was sucnr <dig> .
the adjusted p-values were smaller for tests between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig> than between t= <dig> and t=+ <dig> .
this shows that the transcripts were more differentially expressed between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig> than between t= <dig> and t=+ <dig>  main effects of cortisol released by acth injection on gene expressions are thus observed at t=+ <dig> and t=+ <dig> with a return to baseline levels at t=+ <dig>
genes are divided into clusters corresponding to the kinetics of the response to acth.
full description of the genes including their probe name and localisation is displayed in additional file 1
hac performed on the within-subject deviation matrix with the list of de genes identified  <dig> groups of genes.
figure  <dig> shows that the  <dig> unique de genes allow for an almost perfect classification of the observations with respect to their time of measurement.
for every cluster, fig.
<dig> shows the average evolution of each gene  and the average evolution over the genes in the cluster .
fig.
<dig> hierarchical ascending classification of the  <dig> unique de genes.
a ward method was used with an euclidean distance matrix based on the correlations between genes.
genes are shown in column.
observations are shown in line with one line being a combination pig × time.
colors on the row dendrogram identify the time of measurement.
black: t=0; red: t=+1; green: t=+4; blue: t=+ <dig>  numbers on the column dendrogram identify each cluster
fig.
<dig> average evolution of the genes in each of the cluster identified by the hac on the  <dig> unique de genes.
evolution of each gene is translated so that it is equal to  <dig> at t=0; gray: average evolution of each of the genes in the cluster.
red: average evolution over all genes in the cluster
each cluster was then subjected to a functional analysis .
in each cluster, genes were de  at each time point except for t=+ <dig> in cluster  <dig> .
the first cluster  was characterized by genes increasing with a peak of expression at t= <dig> and stable between t=+ <dig> and t=+ <dig>  the de genes of this cluster could be assembled into a functional network principally involved in neuroimmune functions.
the present analysis reveals novel effects of acth on at least five genes related to immunoregulation .
fkbp <dig>  is a decisive factor for the physiological stress response  <cit>  and has an important role in stress-related phenotypes  <cit> .
it modifies steroid hormone receptor sensitivity  <cit> .
cebpb, dusp <dig>  fkbp <dig> and nfkb1a genes from this cluster are also involved in glucocorticoid receptor signaling.
glucocorticoids exert their classic anti-inflammatory role by acting on nearly all cell types of the immune system.
the ccaat/enhancer binding proteins  are key regulators of cell differentiation and are also involved in the expression and production of inflammatory cytokines  <cit> .
the increase of period  <dig> gene  expression in peripheral blood cells by glucocorticoids was previously reported in humans  <cit> .
physical and inflammatory stressors induce the release of the adrenal glucocorticoid hormone that rapidly alter the expression of per <dig> in peripheral tissues through a gre enhancer present in the gene promotor .
this gene is involved in the circadian rythm, in which the glucocorticoid mechanism plays a predominant role  <cit> .
another de gene ddit <dig>  was described as a surrogate biomarker of the efficiency of glucocorticoid receptor blockade in skeletal muscle  <cit> .
britto and collaborators showed that ddit <dig> expression was low under basal conditions but was highly increased in response to several catabolic stressors, like hypoxia and glucocorticoids  <cit> .
glucocorticoids were shown to up-regulate dusp <dig> in peripheral tissues  <cit>  but constrain the increase of dusp <dig> gene expression in the central components of the hpa axis  <cit> .
in vitro studies have shown that glucocorticoid suppression of some map-kinase dependent cellular processes depends on glucocorticoid mediated up-regulation of dusp <dig> gene expression  <cit> .
the second cluster  was characterized by genes with an increase between t= <dig> and t=+ <dig> and a decrease between t=+ <dig> and t=+ <dig>  this cluster with genes up-regulated at t=+ <dig> is largely related to biological processes such as inflammatory and immune response and genes of which products are located in the plasma membrane.
among these genes, two are particularly interesting.
cd <dig> gene is a component of the innate immune system and has been shown to be sensitive to stress in pigs  <cit> .
megf <dig> gene was shown to be induced by cortisol in human fetal cells in vitro  <cit> .
the third cluster  includes the genes decreasing between t= <dig> and t=+ <dig> and returning to a basal level between t=+ <dig> and t=+ <dig>  no ontology was significantly enriched by genes of this cluster.
it is interesting to underline here the alox <dig> gene  which is a member of the alox family and related to cancer and immune responses.
this gene was also reported as a dexamethasone-responsive gene with nearby glucocorticoid receptor-binding sites  <cit> .
the genes related to the fourth cluster  decrease between t= <dig> and t=+ <dig>  increase between t=+ <dig> and t=+ <dig> and decrease between t=+ <dig> and t=+ <dig>  the fourth cluster corresponds to genes with an overall expression decreasing between t= <dig> and t=+ <dig>  they are significantly linked to biological processes such as protein phosphorylation and kinase activity.
among the genes involved in this cluster arhgap <dig> and arhgap family genes were found to be differentially expressed in macrophages treated with dexamethasone  <cit> .
integration of biological and gene expression data
all de genes were found significantly differentially expressed over time in the mixed model described in section “integration” .
these genes are thus differentially expressed over time even when adjusting for changes in white cells populations.
among them,  <dig> genes had their expression significantly negatively influenced by l/g ratio , meaning that these genes are over-expressed when l/g ratio decreases.
genes with a significant effect of l/g ratio were mainly identified as genes of cluster  <dig>  over-expressed at t=+ <dig>  and cluster  <dig>  over-expressed at t=+ <dig>  and to a lesser extent as genes of cluster  <dig>  under-expressed at t=+ <dig> .
no gene of cluster  <dig> was significantly explained by l/g ratio.
results of the functional analysis of this list of  <dig> genes are shown in additional file  <dig>  biological functions significantly enriched include regulation of apoptotic process, response to lipopolysaccharide, inflammatory and innate immune response, defense response to bacterium and positive regulation of nf-kappab transcription factor activity.
the  <dig> de genes and the biological variables were then subjected to a multilevel pls.
figure  <dig> shows that the first axis of the multilevel pls opposes observations at t=+ <dig> after injection to all others, while the second axis opposes observations at t=+4vs.
all others, similarly as what was already established in multi-level pca of the biological variables in section “overall effect of the injection of exogenous acth on clinical biology variables”.
fig.
<dig> pls regression predicting the biological variables responding to acth from the deg expression a projection of the observations with one point being a combination: time × pig; black: t=0; red: t=+1; green: t=+4; blue: t=+ <dig>  b projection of the variables; blue: biological variables; red: gene expressions;  <dig> genes/components were kept using a sparse approach
on the first axis, cortisol and ffa levels are strongly positively correlated with the expressions of cebpb, rgs <dig>  rhob, per <dig>  fkbp <dig>  cebpd, ddit <dig>  cpt1a and dusp <dig>  all these genes belong to the first cluster identified earlier and are linked to molecular functions such as protein binding and transcription regulation.
the second axis of the multilevel pls is characterized by the opposition between the proportion of lymphocytes and monocytes vs. the proportion of granulocytes.
this axis is positively correlated with sucnr <dig>  slco2b <dig>  fbp <dig> and loc <dig>  these genes belong to the third cluster and are related to glycolysis and glycogenesis.
sucnr <dig>  is decreased at t=+ <dig> and increased at t=+ <dig>  succinate has a wide range of metabolic actions and regulates the functions of macrophages  <cit> .
the axis is negatively correlated with cd <dig>  clc4d, chit <dig>  megf <dig> and c2h19orf <dig>  these genes belong to the second cluster which is linked to molecular functions such as inflammatory response, but their relationships with cortisol or stress are not yet clearly established.
eight genes  are functionally connected to each other by nr3c <dig>  the nr3c <dig>  is the glucocorticoid receptor, which can function both as a transcription factor that binds to glucocorticoid response elements in the promoters of glucocorticoid responsive genes, and as a regulator of other transcription factors.
functional consequences of glucocorticoid receptor polymorphisms were reported in pigs .
mutations in nr3c <dig> have been previously demonstrated to be associated with generalized glucocorticoid resistance  <cit> .
it is interesting to highlight the de genes that encode transcription factors.
they play a crucial role in regulating gene expression and are fit to regulate diverse cellular processes by interacting with other proteins.
most of them have not yet been described as important in transcription networks involved in stress responses.
if the genes are co-expressed it is highly probable that they are co-regulated.
this knowledge can provide new patterns of biomarkers of the individual sensitivity to cortisol that is our field of interest in this study.
fig.
<dig> gene network related to glucocorticoid response in whole blood transcriptome  <dig> h after acth injection this network corresponds to the genes up-regulated  <dig> h after acth injection .
it combines bibliographic  and regulatory relationships , green lines) proposed by ingenuity software.
cortisol has the highest plasma level at t=+ <dig> and acts through the glucocorticoid receptor nr3c1
our results are in accordance with several studies on the effects of glucocorticoid hormones on peripheral blood cells.
numerous genes related to cluster  <dig> and shown as acth responsive were found differentially expressed in stress-related investigations.
five genes found in our study  were proposed as markers of differential glucocorticoid sensitivity  <cit> .
nfkbia, dusp <dig>  cebpd, fkbp <dig> genes were also found to be associated with up- and down-regulated clusters in response to continuous  <dig> h cortisol infusion  <cit> .
ponsuksili and collaborators  <cit>  describe nfkbia, cebpb and cebpd as genes of which hepatic expression levels are correlated with plasma cortisol concentrations.
up-regulation of per <dig> gene upon gr activation was confirmed by genome-wide study of glucocorticoid receptor binding sites in neuronal pc <dig> cells  <cit> .
however, ddit <dig> was shown to be down-regulated by gr activation rather than up-regulated in this analysis.
while looking for genes of which expressions at t= <dig> t=+ <dig> t=+ <dig> or t=+ <dig> were significantly correlated with the level of cortisol at t=+ <dig>  only two genes were identified: trmt2a , a gene involved in the methylation of trna, and loc <dig> , a gene of which function has not been identified yet.
there is a negative relationship between the expression of cortisol at t=+ <dig> and the expression of these two genes at t= <dig> .
this implies that when their baseline expression is higher, the intensity of the cortisol response to acth decreases.
CONCLUSIONS
the present work shows the interest of transcriptomic data analysis at multiple levels.
in other studies, genetic markers found through an analysis of transcription factor binding sites of differentially expressed genes in peripheral blood cells have been proposed in humans to identify the chronic stress related to psychopathological conditions  <cit> .
in farm animals, this approach was used in horses  <cit> .
these studies show chronic stress-related changes in the balance between the expression of stress-related genes regulated by glucocorticoids and those regulated by inflammation-related factors.
furthermore, recent data in humans show that the immune system function can also be assessed through blood transcriptomics in health and disease  <cit> .
in the present study, we identified  <dig> genes differentially expressed in peripheral blood cells of pigs in response to acth at different times after injection.
it therefore supplies biological markers of hpa axis activation at the gene expression level, and the knowledge on functional gene clusters will help to elucidate the biological processes involved.
moreover, these genes might be candidates for a better understanding of the molecular mechanisms related to stress responses.
thus, blood transcriptome analysis appears as a promising avenue to develop multidimensional biological markers related to robustness.
these markers should be used in the study of the genetic mechanisms of adaptation in farm animals that will help to deliver genetic strategies to animal breeders in order to balance production objectives and robustness of animals as well as their welfare  <cit> .
availability of supporting data
the data sets supporting the results of this article are available in gene expression omnibus  through the accession number gse71207).
additional files
additional file  <dig> ∙ gene name: name of the gene
∙ probe name, agilent gpl16524: probe names used on the microarray agilent gpl16524
∙ l/g coefficient: l/g ratio coefficient estimate
∙ adjusted pval : adjusted p-value associated with the l/g ratio coefficient
∙ adj.pval: adjusted p-value of the test at the time measurement where the most significant duplicate of the gene is de
∙ time point: time measurement where the gene is de
∙ expression: whether the deg is up or down-regulated
∙ cluster: cluster in which the gene is classified by hac
∙ gene description: informations on the gene’s molecular function
∙ location: chromosomal location of the gene.
additional file  <dig>
distribution of the rank of the significant adjusted
p
-values in the tests for de transcripts between
t
=  <dig> and
t
= +  <dig>
t
=  <dig> and
t
= +  <dig> and
t
=  <dig> and
t
= +  <dig> ‘.pdf’ file.
p-values are smaller at t=+ <dig> and t=+ <dig> than at t=+ <dig> implying that the transcripts were overall more differentially expressed between t= <dig> and t=+ <dig> and between t= <dig> and t=+ <dig> than between t= <dig> and t=+ <dig>
additional file  <dig>
complete list of enriched go , molecular function  and cellular component  for each of the cluster identified with the hierarchical ascending clustering ‘.xls’ file.
features the go items, the corresponding functions, the class of ontology, the number of genes in the input list  and in the reference list , the raw and the adjusted fisher’s exact test p-value and the list of genes.
additional file  <dig>
complete list of enriched go , molecular function  and cellular component  for  <dig> genes for which l/g ratio had a significant effect ‘.xls’ file.
features the go items, the corresponding functions, the class of ontology, the number of genes in the input list  and in the reference list , the raw and the  <dig> adjusted fisher’s exact test p-value and the list of genes.
abbreviations
acthadrenocorticotropic hormone
bhbenjamini-hochberg
bwbirth weight
cbgcorticosteroid binding globulin
dedifferentially expressed
fdrfalse discovery rate
grglucocorticoid receptor
greglucocorticoid response element
hachierarchical ascendant classification
hcthematocrit
hgbhemoglobin
hpahypothalamic-pituitary-adrenocortical
l/ glymphocytes/ granulocytes
pcaprinicpal component analysis
plspartial least square regression
rcred cells
semstandar error of the mean
wwweaning weight
competing interests
the authors declare that they have no competing interests.
authors’ contributions
vs and nvv developed and performed all the statistics.
pm and et developed and undertook the experimental design and ensured the biological interpretation.
lg, et and ll ensured the experimental transcriptomic analysis.
yl provided the transcriptomic data set vs, et and ll performed the transcriptomic related biology interpration.
lg performed the sampling and the data management.
cl performed the genetics and the choice of animals.
yb was responsible of the animal management.
vs, nvv, et, pm and ll wrote the paper.
pm is the project manager.
all authors read and approved the final manuscript.

 farm animals have been highly selected for favorable production traits such as lean growth rate, feed efficiency, and prolificacy in pigs. in parallel their robustness has tended to decrease, as shown by the sensitivity to diseases, locomotor weakness or behavioral problems. the hypothalamic-pituitary-adrenocortical axis plays a major role in physiological homeostasis including metabolism, brain function and behavior, the immune system and inflammatory processes. together with the autonomic nervous system, it is also involved in stress and adaptive responses to environmental challenges. on the basis of available literature, we have hypothesized that the hpa axis could contribute to the trade-off between production and robustness traits, and that genetic variation in hpa axis activity could be used to select more robust animals. hpa axis activity shows a large variation among individuals and genetic influences are well documented. for example, in pigs, the sensitivity of the adrenal cortex to adrenocorticotropic hormone and the production of corticosteroid binding globulin, the carrier of cortisol in blood are the two main mechanisms responsible for genetic differences in circulating cortisol levels. in a previous paper, hazard et al. have studied at the gene expression level the molecular mechanisms of genetic differences in the adrenal gland response to acth. however little is known about the individual variation in the biological activity of cortisol, the main glucocorticoid hormone, and the genetic mechanisms involved. corticosteroid hormones exert their biological actions via two intracellular receptors that, upon activation by their ligand, influence the expression of a large number of genes in a wide range of cell types. in pigs, glucocorticoid receptor polymorphisms have been described with their functional consequences. the present experiment was designed to explore in pigs individual variation in the biological response to cortisol in order to identify possible biomarkers of this sensitivity. to that end, juvenile pigs were submitted to an acth challenge to release cortisol and the kinetics of the response was measured with biological variables and with gene expression analysis in blood cells. taken together with biological information, this approach will serve as an important step to understand hpa axis regulation and will identify key genes involved in signaling pathways relevant to stress responses. the final goal of this work is to develop a strategy for further genetic studies in order to overcome the unfavourable consequences of stress in farm animals. animals and methods animals, treatment and blood sampling all animal use was performed under european union and french legislation. the protocol and procedures were approved by the local ethics committee. experimental animals were piglets randomly selected from litters of purebred large white pigs and produced in successive batches raised weeks apart. they were weaned at weeks and animals from 2– litters were mixed at weaning in each post-weaning pen. experimental animals were not isolated from their littermates. at weeks, each animal was injected in the neck muscles with μg of synthetic 1– acth diluted in ml of % saline. injections occurred from 10:00–11: am to avoid nycthemeral variations. blood samples were collected before the injection and, and hours later. at each time, animals were slightly restrained on their back in such a way that the effect on their stress level can be regarded as insignificant. two blood samples were then taken by puncture of a jugular vein in vacutainer Ⓡ tubes with g needles. the whole handling procedure lasted less than sec. one ml tube with lithium heparin was used for chemical biology. after centrifugation, plasma aliquots were frozen at – °c until analysis. one ml tube with edta was used for blood cell count and an aliquot was mixed with the same volume of dl buffer, frozen at – °c for h and then at – °c until analysis for gene expression. biological analyses cortisol was measured by direct automated immunoassay. glucose and free fatty acid, were measured by colorimetry with an abx pentra clinical chemistry analyzer from horiba medical. blood cell counts were measured with a ms-9- hematology analyzer from melet schloesing laboratories, calibrated for pig blood by the manufacturer. blood cell count variables included: white cells count, proportion of lymphocytes, monocytes and granulocytes, red cells count, percentage of hematocrit, concentration of hemoglobin, red cells width and volume, concentration of platelets and platelets width and volume. the biological variables contained thus variables measured on the pigs. in addition, birth and weaning weights were also measured for each pig. outlying observations were visually identified and treated as missing data. missing data were imputed using a k-nearest neighbour imputation. to ensure normality, cortisol, platelets and white cells count were log transformed and ffa was transformed using the square-root. batch effects were removed by aligning the within-batch medians for all measurements. rna extraction and whole blood gene expression analysis a total rna isolation and purification was done according to the manufacturer’s instructions using the nucleospin rna blood kit followed by dnase treatment. the quality of each rna sample was checked through the bioanalyser agilent and low-quality rna preparations were discarded. microarray description a porcine microarray gpl was used to hybridize the rna samples as previously described. this microarray contained spots. among them, were negative controls and were used for aligning. one probe was duplicated twice on each array. thus the microarray contained unique porcine probes. after quality control, quantile normalization and filtering, transcripts were found to be expressed in blood in our experimental conditions. hybridization protocol blood samples from female pigs from only batches were used. each of the arrays used contained microarrays which corresponded to the observations of two individuals, each from one batch. this design secured the kinetics of the response for each individual and prevented confounding effects between batch and array. after quality control and filtering, probes were kept and log transformed. technical biases were handled by aligning the within-array medians for all genes and by a quantile normalization within animal. no missing data were reported. normalized data were submitted to geo/ncbi with the gse accession number. statistical analyses all analyses were performed with the r software, version.. multilevel approach a multilevel approach was used to investigate the relationships between the repeated measurements while taking advantage of multivariate approaches. the multilevel approach, as described by liquet et al., uses a split-up variation coming from the mixed-model framework. let x=i= …,n,t∈{ + + +24},k= …,p be the observation matrix on n animals with times of measurements. x can be split up as: x=x..⏟offset term+xb⏟between-animal deviation+xw⏟within-animal deviation the matrix x.. represents the offset term defined as 1nx..t where 1n is a vector of length n containing only ones and x..t=. xb is the between-animal matrix of size defined by concatenating 14xbit for each animal into xb with xbit=xi.1−x.. …,xi.p−x..pxi.k=14∑t∈{ + + +24}xitk. xw=x−xi. is the within-animal deviation matrix of size with xi. the matrix defined by concatenating the matrices 14xi.t for every animal into xi., with xi.t=xi. …,xi.p. by splitting the different parts of the variation in the data while taking into account the repeated measurements on the subjects, the multilevel step allows us to study the effect of different conditions within a subject separately from the variation between subjects. this method is especially relevant when a high between-subject variability is observed in repeated data: multivariate approaches can then be performed on xw to highlight the most relevant correlations between variables in the dataset, independently from individual variations. statistical analysis of plasma metabolites and hormone first, all variables were subjected to a one-way anova with repeated measures. p-values were adjusted using a benjamini-hochberg approach in order to control the false discovery rate. variables with an adjusted p-value were then subjected to paired t-tests to assess the difference between t= and t=+ between t= and t=+ and between t= and t=+ the full list of p-values was adjusted using a bh approach. in addition, the influence of sex on the biological variables was tested using a two-way anova with repeated measures including sex as a variable. p-values were adjusted using a bh approach. then, variability between individuals before the acth injection was studied by performing a pca on the observations at t= for variables responding to acth and birth and weaning weights. the overall effect of acth over time was investigated with a multilevel pca as previously described. cortisol levels at t=+ is the most relevant measure to assess the sensitivity of the adrenals to acth. hence, correlations between biological variables at t∈{ + + +24} and the level of cortisol at t=+ were investigated using paired t-tests. p-values were adjusted using a bh approach. statistical analysis of the transcriptome all transcripts were subjected to paired t-tests to assess the difference in expression between t= and t=+ between t= and t=+ and between t= and t=+ the full list of p-values was adjusted using a bonferroni approach. as the bonferroni approach exerts a more stringent control than the bh approach, it was used to obtain a narrowed list of the most significantly differentially expressed transcripts. transcripts with at least one adjusted p-value < among the three tests were kept. correlations between the expression levels of different transcripts of the same gene were investigated to highlight genes for which at least 3/ of the duplicates were significantly de and had a correlation of at least between them. the most significant transcripts per annotated gene were kept and the multilevel approach was used to extract the within-subject deviation matrix for further analysis. a hierarchical ascending classification was performed using the ward method with an euclidean distance matrix based on the correlations between genes. this allowed for the characterization of clusters of genes having the same evolution over time. significance of the difference between time measurements was assessed using paired t-tests at the average gene level between t= and t=+ between t= and t=+ and between t= and t=+ p-values were adjusted within the clusters using a bh approach. integration relations between the main stress variable, the cortisol, the biological variables and the gene expressions were studied using different approaches. as for biological data, correlations between gene expression at t∈{ + + +24} and the level of cortisol at t=+ were investigated. more precisely, the within-subject components of the transcriptomic data and of the cortisol level was extracted using the method described in and section “multilevel approach”. then, a paired t-test was performed between the within-subject expression at t∈{ + + +24} and the within-subject level of cortisol at t=+ the full list of p-values was globally adjusted using a bh approach. in order to assess the significativity of change in expression after removing the contribution of cell population changes, a linear mixed model was fitted for every deg xit=β0+β t+β2it+ui+εit where xit is the expression of the deg for the animal number i and the time step t, it is the lymphocytes/granulocytes ratio for the same experiment and ui is the individual random effect. both time step and it were supposed to have fixed effects on gene expression. significance of the time effect in this model was checked by testing β t=β for t∈{+ + +24} and correcting p-values by time point for multiple tests with a bh approach. the effect of changes in white cell populations was also assessed by testing β2= and a correction for multiple tests was applied using a bh approach. a multilevel pls, i.e., a pls performed on the within-subject components of the biological and transcriptomic data, as described in and section “multilevel approach”, was used to investigate the overall relationships between biological and transcriptomic data. a sparse version of the pls as described in was used to select a small subset of variables to explain every axis. sequence annotation each cdna sequence was compared to refseq_rna mammalian database using the ncbi blastn program. resulting hits were sorted out according to their closeness to the pig genome, their coverage and sequence identity. the selected cdna sequences were submitted to the hugo gene nomenclature committee, using their refseq ids. then, hugo gene symbols or official gene symbol were used as gene names. functional enrichment and pathway analysis functional enrichment was performed on each list of clustered genes identified by hac and on the list of genes of which expression is significatively explained by cell population changes in blood. functional annotation of genes was provided by the biomart database. to set the statistical enrichment of a particular biological function, a fisher’s exact test was performed, using the list of genes expressed on the microarray as the reference list of genes. resulting p-values were adjusted for multiple tests using a bh approach. a minimum of genes per gene ontology and a fdr< were necessary to consider a biological function to be enriched. a pathway is an interconnected arrangement of processes, representing the functional roles of genes in the genome. functional integration of gene expression, i.e., identification of gene networks, was performed using the ‘gene ontology’ database amigo. the significantly up- or down-regulated genes could be assembled into networks using ‘ingenuity pathway analysis’ under licence. this application provides computational algorithms to identify the enriched biological pathways, functions and biological mechanisms of selected genes and proposes also enriched regulators as transcription factors. a regulatory network could be constructed with the information provided by the option ‘upstream regulator’. this option proposes a list of regulators known to have a significant effect on some of the targeted genes in the input list. ingenuity also provides computational algorithms to identify and to dynamically generate significant biological networks. networks are ranked by a score that takes into account the number of focus genes and the size of the networks. this score ) indicates the probability for genes to be associated in the same network by chance. a higher score means a smaller probability for genes to be observed in the same network by chance. we chose networks displaying direct relationships between genes. path designer was used to improve the readability of the networks. nodes added by ingenuity were discarded when they were not necessary to connect our genes of interest and the resulting network was merged with the regulatory network.  plasma cortisol, metabolites and hematologyfig. mean evolution of the biological variables overtime.: time measurement at which the expression of the variable is significantly different from the expression at t= vertical bars correspond to + and – sem at each time point f value and fdr are for the test of the global time effect on each variable. nr non relevant since the measure is the same at all time steps as expected, acth induced a strong cortisol response peaking h after injection. plasma cortisol levels at t=+ were -fold higher than and highly correlated with basal levels. they were significantly lower than baseline values at t=+ and at t=+, as can be expected from the feedback of cortisol on its own secretion. plasma glucose levels showed a slight increase at t=+ after acth injection, but ffa levels showed a large increase at t=+, with a strong variation across animals, and were almost back to basal levels at t=+ the values measured at t=+ were correlated with basal values. cortisol induces a weak increase in circulating glucose and also potentiates the effect of other counter-regulatory hormones and increases ffa levels via acute lipolysis. the data obtained here for clinical hematology measures are in the range of published values in pigs. these values show large variations with age and breed among other sources. although the total number of leucocytes was only marginally influenced by acth, large changes in leucocyte subpopulations were observed with an increase of the proportion of granulocytes and a decrease of the proportion of lymphocytes and monocytes at t=+ and t=+ these effects are consistent with the literature and result from the redistribution of leucocytes between blood and tissues. red blood cell related variables were decreased after injection of acth and remained low at t=+ platelets were not influenced by injection of acth. sex did not influence any of the variables. between-subject variability at t=0 results of the pca at t= on variables responding to acth and birth and weaning weights are shown in fig. as red cells count, hematocrit and hemoglobin are redundant variables, decision was made to keep only rc for the pca. no outliers are identified on the projection of the individuals on the two first dimensions of the pca. on this pca, sex was not found to be related to the main variability in the dataset. the first dimension shows an opposition between the proportion of granulocytes and the proportion of lymphocytes. the second axis shows an opposition between cortisol and birth and weaning weights with a strong opposition between these variables on the whole plan formed by the first and second axis. the other variables were not correlated with either of the first two dimensions of the pca. fig. pca on the biological variables identified as responding to acth, the birth and weaning weights at t= colors symbolize the sex: black = female; red = male a projection of the individuals on dimensions and 2; b projection of the variables on dimensions and 2; bw = birth weight; ww = weaning weight; lympho = lymphocyte ratio; mono = monocytes ratio; granulo = granulocyte ratio; rc = red cell count; gluc = glucose; ffa = free fatty acids overall effect of the injection of exogenous acth on clinical biology variables extraction of the within-subject data matrix prior to the application of a pca analysis allows for the separation of the observations according to their time of measurement. the first component of the multi-level pca opposes the observations at t= and t=+ to the observations at t=+, this time step corresponding to the peak of cortisol. the second component opposes the observations at t=+ to the observations at t=+. the representation of the variables shows that the first axis is mainly driven by an opposition between the proportion of granulocytes, ffa, cortisol and red cell count, on one side, and lymphocytes and monocytes, on the other side. the second axis shows an opposition between glucose and cortisol, ffa and red count. fig. multilevel pca on the biological variables responding to acth. colors symbolize the time of measurement; black: t=0; red: t=+1; green: t=+4; blue: t=+24; a projection of the individuals on dimensions 1–2; b projection of the variables on dimensions 1–2; lympho = lymphocyte ratio; mono = monocyte ratio; granulo = granulocyte ratio; rc = red cell counts; gluc = glucose; ffa = free fatty acids specific links to the level of cortisol at t=+1 correlations between cortisol at t=+ and other variables at t= allows for the identification of variables which baseline levels may be directly or indirectly linked to the intensity of the cortisol level in response to acth, a measure of individual variation in hpa axis activity. correlations are shown in table only glucose and ffa levels at t= were significantly positively correlated with the level of cortisol at t=+. correlations between cortisol at t=+ and other variables at t=+ t=+ or t=+ allows for the identification of variables which are directly linked to the level of cortisol when it reaches its peak during the stress response. ffa at t=+ were positively correlated with cortisol at t=+ and glucose at t=+ t=+ and t=+ was negatively correlated with cortisol at t=+. no other variable was found to be significantly linked to the intensity of the cortisol level in response to acth. – – – se: standard error of the correlation coefficient; in bold: significantly ≠ differentially expressed genes we used a comprehensive gene expression profiling by means of microarray analysis to identify clusters of genes differentially expressed in peripheral blood cells, taking into consideration the kinetics of the response with time points. differential analysis revealed de transcripts matching unique genes. among them, genes were differentially expressed at t=+, were differentially expressed at t=+ and were differentially expressed at t=+. the only gene de at both t=+ and t=+ was sucnr. the adjusted p-values were smaller for tests between t= and t=+ and between t= and t=+ than between t= and t=+. this shows that the transcripts were more differentially expressed between t= and t=+ and between t= and t=+ than between t= and t=+ main effects of cortisol released by acth injection on gene expressions are thus observed at t=+ and t=+ with a return to baseline levels at t=+ genes are divided into clusters corresponding to the kinetics of the response to acth. full description of the genes including their probe name and localisation is displayed in additional file 1 hac performed on the within-subject deviation matrix with the list of de genes identified groups of genes. figure shows that the unique de genes allow for an almost perfect classification of the observations with respect to their time of measurement. for every cluster, fig. shows the average evolution of each gene and the average evolution over the genes in the cluster. fig. hierarchical ascending classification of the unique de genes. a ward method was used with an euclidean distance matrix based on the correlations between genes. genes are shown in column. observations are shown in line with one line being a combination pig × time. colors on the row dendrogram identify the time of measurement. black: t=0; red: t=+1; green: t=+4; blue: t=+ numbers on the column dendrogram identify each cluster fig. average evolution of the genes in each of the cluster identified by the hac on the unique de genes. evolution of each gene is translated so that it is equal to at t=0; gray: average evolution of each of the genes in the cluster. red: average evolution over all genes in the cluster each cluster was then subjected to a functional analysis. in each cluster, genes were de at each time point except for t=+ in cluster. the first cluster was characterized by genes increasing with a peak of expression at t= and stable between t=+ and t=+ the de genes of this cluster could be assembled into a functional network principally involved in neuroimmune functions. the present analysis reveals novel effects of acth on at least five genes related to immunoregulation. fkbp is a decisive factor for the physiological stress response and has an important role in stress-related phenotypes. it modifies steroid hormone receptor sensitivity. cebpb, dusp fkbp and nfkb1a genes from this cluster are also involved in glucocorticoid receptor signaling. glucocorticoids exert their classic anti-inflammatory role by acting on nearly all cell types of the immune system. the ccaat/enhancer binding proteins are key regulators of cell differentiation and are also involved in the expression and production of inflammatory cytokines. the increase of period gene expression in peripheral blood cells by glucocorticoids was previously reported in humans. physical and inflammatory stressors induce the release of the adrenal glucocorticoid hormone that rapidly alter the expression of per in peripheral tissues through a gre enhancer present in the gene promotor. this gene is involved in the circadian rythm, in which the glucocorticoid mechanism plays a predominant role. another de gene ddit was described as a surrogate biomarker of the efficiency of glucocorticoid receptor blockade in skeletal muscle. britto and collaborators showed that ddit expression was low under basal conditions but was highly increased in response to several catabolic stressors, like hypoxia and glucocorticoids. glucocorticoids were shown to up-regulate dusp in peripheral tissues but constrain the increase of dusp gene expression in the central components of the hpa axis. in vitro studies have shown that glucocorticoid suppression of some map-kinase dependent cellular processes depends on glucocorticoid mediated up-regulation of dusp gene expression. the second cluster was characterized by genes with an increase between t= and t=+ and a decrease between t=+ and t=+ this cluster with genes up-regulated at t=+ is largely related to biological processes such as inflammatory and immune response and genes of which products are located in the plasma membrane. among these genes, two are particularly interesting. cd gene is a component of the innate immune system and has been shown to be sensitive to stress in pigs. megf gene was shown to be induced by cortisol in human fetal cells in vitro. the third cluster includes the genes decreasing between t= and t=+ and returning to a basal level between t=+ and t=+ no ontology was significantly enriched by genes of this cluster. it is interesting to underline here the alox gene which is a member of the alox family and related to cancer and immune responses. this gene was also reported as a dexamethasone-responsive gene with nearby glucocorticoid receptor-binding sites. the genes related to the fourth cluster decrease between t= and t=+ increase between t=+ and t=+ and decrease between t=+ and t=+ the fourth cluster corresponds to genes with an overall expression decreasing between t= and t=+ they are significantly linked to biological processes such as protein phosphorylation and kinase activity. among the genes involved in this cluster arhgap and arhgap family genes were found to be differentially expressed in macrophages treated with dexamethasone. integration of biological and gene expression data all de genes were found significantly differentially expressed over time in the mixed model described in section “integration”. these genes are thus differentially expressed over time even when adjusting for changes in white cells populations. among them, genes had their expression significantly negatively influenced by l/g ratio, meaning that these genes are over-expressed when l/g ratio decreases. genes with a significant effect of l/g ratio were mainly identified as genes of cluster over-expressed at t=+ and cluster over-expressed at t=+ and to a lesser extent as genes of cluster under-expressed at t=+. no gene of cluster was significantly explained by l/g ratio. results of the functional analysis of this list of genes are shown in additional file biological functions significantly enriched include regulation of apoptotic process, response to lipopolysaccharide, inflammatory and innate immune response, defense response to bacterium and positive regulation of nf-kappab transcription factor activity. the de genes and the biological variables were then subjected to a multilevel pls. figure shows that the first axis of the multilevel pls opposes observations at t=+ after injection to all others, while the second axis opposes observations at t=+4vs. all others, similarly as what was already established in multi-level pca of the biological variables in section “overall effect of the injection of exogenous acth on clinical biology variables”. fig. pls regression predicting the biological variables responding to acth from the deg expression a projection of the observations with one point being a combination: time × pig; black: t=0; red: t=+1; green: t=+4; blue: t=+ b projection of the variables; blue: biological variables; red: gene expressions; genes/components were kept using a sparse approach on the first axis, cortisol and ffa levels are strongly positively correlated with the expressions of cebpb, rgs rhob, per fkbp cebpd, ddit cpt1a and dusp all these genes belong to the first cluster identified earlier and are linked to molecular functions such as protein binding and transcription regulation. the second axis of the multilevel pls is characterized by the opposition between the proportion of lymphocytes and monocytes vs. the proportion of granulocytes. this axis is positively correlated with sucnr slco2b fbp and loc these genes belong to the third cluster and are related to glycolysis and glycogenesis. sucnr is decreased at t=+ and increased at t=+ succinate has a wide range of metabolic actions and regulates the functions of macrophages. the axis is negatively correlated with cd clc4d, chit megf and c2h19orf these genes belong to the second cluster which is linked to molecular functions such as inflammatory response, but their relationships with cortisol or stress are not yet clearly established. eight genes are functionally connected to each other by nr3c the nr3c is the glucocorticoid receptor, which can function both as a transcription factor that binds to glucocorticoid response elements in the promoters of glucocorticoid responsive genes, and as a regulator of other transcription factors. functional consequences of glucocorticoid receptor polymorphisms were reported in pigs. mutations in nr3c have been previously demonstrated to be associated with generalized glucocorticoid resistance. it is interesting to highlight the de genes that encode transcription factors. they play a crucial role in regulating gene expression and are fit to regulate diverse cellular processes by interacting with other proteins. most of them have not yet been described as important in transcription networks involved in stress responses. if the genes are co-expressed it is highly probable that they are co-regulated. this knowledge can provide new patterns of biomarkers of the individual sensitivity to cortisol that is our field of interest in this study. fig. gene network related to glucocorticoid response in whole blood transcriptome h after acth injection this network corresponds to the genes up-regulated h after acth injection. it combines bibliographic and regulatory relationships, green lines) proposed by ingenuity software. cortisol has the highest plasma level at t=+ and acts through the glucocorticoid receptor nr3c1 our results are in accordance with several studies on the effects of glucocorticoid hormones on peripheral blood cells. numerous genes related to cluster and shown as acth responsive were found differentially expressed in stress-related investigations. five genes found in our study were proposed as markers of differential glucocorticoid sensitivity. nfkbia, dusp cebpd, fkbp genes were also found to be associated with up- and down-regulated clusters in response to continuous h cortisol infusion. ponsuksili and collaborators describe nfkbia, cebpb and cebpd as genes of which hepatic expression levels are correlated with plasma cortisol concentrations. up-regulation of per gene upon gr activation was confirmed by genome-wide study of glucocorticoid receptor binding sites in neuronal pc cells. however, ddit was shown to be down-regulated by gr activation rather than up-regulated in this analysis. while looking for genes of which expressions at t= t=+ t=+ or t=+ were significantly correlated with the level of cortisol at t=+ only two genes were identified: trmt2a, a gene involved in the methylation of trna, and loc, a gene of which function has not been identified yet. there is a negative relationship between the expression of cortisol at t=+ and the expression of these two genes at t=. this implies that when their baseline expression is higher, the intensity of the cortisol response to acth decreases.  the present work shows the interest of transcriptomic data analysis at multiple levels. in other studies, genetic markers found through an analysis of transcription factor binding sites of differentially expressed genes in peripheral blood cells have been proposed in humans to identify the chronic stress related to psychopathological conditions. in farm animals, this approach was used in horses. these studies show chronic stress-related changes in the balance between the expression of stress-related genes regulated by glucocorticoids and those regulated by inflammation-related factors. furthermore, recent data in humans show that the immune system function can also be assessed through blood transcriptomics in health and disease. in the present study, we identified genes differentially expressed in peripheral blood cells of pigs in response to acth at different times after injection. it therefore supplies biological markers of hpa axis activation at the gene expression level, and the knowledge on functional gene clusters will help to elucidate the biological processes involved. moreover, these genes might be candidates for a better understanding of the molecular mechanisms related to stress responses. thus, blood transcriptome analysis appears as a promising avenue to develop multidimensional biological markers related to robustness. these markers should be used in the study of the genetic mechanisms of adaptation in farm animals that will help to deliver genetic strategies to animal breeders in order to balance production objectives and robustness of animals as well as their welfare. availability of supporting data the data sets supporting the results of this article are available in gene expression omnibus through the accession number gse71207). additional files additional file ∙ gene name: name of the gene ∙ probe name, agilent gpl16524: probe names used on the microarray agilent gpl16524 ∙ l/g coefficient: l/g ratio coefficient estimate ∙ adjusted pval: adjusted p-value associated with the l/g ratio coefficient ∙ adj.pval: adjusted p-value of the test at the time measurement where the most significant duplicate of the gene is de ∙ time point: time measurement where the gene is de ∙ expression: whether the deg is up or down-regulated ∙ cluster: cluster in which the gene is classified by hac ∙ gene description: informations on the gene’s molecular function ∙ location: chromosomal location of the gene. additional file distribution of the rank of the significant adjusted p -values in the tests for de transcripts between t = and t = + t = and t = + and t = and t = + ‘.pdf’ file. p-values are smaller at t=+ and t=+ than at t=+ implying that the transcripts were overall more differentially expressed between t= and t=+ and between t= and t=+ than between t= and t=+ additional file complete list of enriched go, molecular function and cellular component for each of the cluster identified with the hierarchical ascending clustering ‘.xls’ file. features the go items, the corresponding functions, the class of ontology, the number of genes in the input list and in the reference list, the raw and the adjusted fisher’s exact test p-value and the list of genes. additional file complete list of enriched go, molecular function and cellular component for genes for which l/g ratio had a significant effect ‘.xls’ file. features the go items, the corresponding functions, the class of ontology, the number of genes in the input list and in the reference list, the raw and the adjusted fisher’s exact test p-value and the list of genes. abbreviations acthadrenocorticotropic hormone bhbenjamini-hochberg bwbirth weight cbgcorticosteroid binding globulin dedifferentially expressed fdrfalse discovery rate grglucocorticoid receptor greglucocorticoid response element hachierarchical ascendant classification hcthematocrit hgbhemoglobin hpahypothalamic-pituitary-adrenocortical l/ glymphocytes/ granulocytes pcaprinicpal component analysis plspartial least square regression rcred cells semstandar error of the mean wwweaning weight competing interests the authors declare that they have no competing interests. authors’ contributions vs and nvv developed and performed all the statistics. pm and et developed and undertook the experimental design and ensured the biological interpretation. lg, et and ll ensured the experimental transcriptomic analysis. yl provided the transcriptomic data set vs, et and ll performed the transcriptomic related biology interpration. lg performed the sampling and the data management. cl performed the genetics and the choice of animals. yb was responsible of the animal management. vs, nvv, et, pm and ll wrote the paper. pm is the project manager. all authors read and approved the final manuscript. 
hpa axis plays a major role in physiological homeostasis. microarray analysis identified unique genes of which expression responded to the injection of acth. a multilevel statistical analysis was used to take into account the longitudinal aspect of the data. these genes might be candidates for a better understanding of the molecular mechanisms of the stress response. basal level of cortisol was correlated with birth and weaning weights. in farm animals in general and specifically in pigs, breeding strategies have highly favored production traits such as lean growth rate, feed efficiency and prolificacy at the cost of robustness. in particular, ddit dusp fkbp il7r, nfkbia, per rgs and rhob were shown to be connected to each other by the glucocorticoid receptor nr3c most of the differentially expressed genes that encode transcription factors have not been described yet as being important in transcription networks involved in stress response. we used acth injections to trigger production of cortisol in juvenile large white pigs from litters and the kinetics of the response was measured with biological variables and whole blood gene expression at time points. cortisol level reached its peak h after acth injection. these genes were classified into clusters with distinctive kinetics in response to acth injection. we identified genes as biological markers of hpa axis activation at the gene expression level. their co-expression may mean co-regulation and they could thus provide new patterns of biomarkers of the individual sensitivity to cortisol. on the hypothesis that the hpa axis could contribute to the trade-off between robustness and production traits, we have designed this experiment to explore individual variation in the biological response to the main stress hormone, cortisol, in pigs.
 hpa axis plays a major role in physiological homeostasis. it is also involved in stress and adaptive response to the environment. in farm animals in general and specifically in pigs, breeding strategies have highly favored production traits such as lean growth rate, feed efficiency and prolificacy at the cost of robustness. on the hypothesis that the hpa axis could contribute to the trade-off between robustness and production traits, we have designed this experiment to explore individual variation in the biological response to the main stress hormone, cortisol, in pigs. we used acth injections to trigger production of cortisol in juvenile large white pigs from litters and the kinetics of the response was measured with biological variables and whole blood gene expression at time points. a multilevel statistical analysis was used to take into account the longitudinal aspect of the data.  cortisol level reached its peak h after acth injection. white blood cell composition was modified with a decrease of lymphocytes and monocytes and an increase of granulocytes. basal level of cortisol was correlated with birth and weaning weights. microarray analysis identified unique genes of which expression responded to the injection of acth. these genes were classified into clusters with distinctive kinetics in response to acth injection. the first cluster identified genes strongly correlated to cortisol and previously reported as being regulated by glucocorticoids. in particular, ddit dusp fkbp il7r, nfkbia, per rgs and rhob were shown to be connected to each other by the glucocorticoid receptor nr3c most of the differentially expressed genes that encode transcription factors have not been described yet as being important in transcription networks involved in stress response. their co-expression may mean co-regulation and they could thus provide new patterns of biomarkers of the individual sensitivity to cortisol.  we identified genes as biological markers of hpa axis activation at the gene expression level. these genes might be candidates for a better understanding of the molecular mechanisms of the stress response. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords stresshypothalamic-pituitary-adrenal axiscortisoltime-coursesystems biologymicroarraypigissue-copyright-statement© the author 2015 