Time course of the response to ACTH in pig: biological and transcriptomic study

Plasma cortisol, metabolites and hematology

Table 1 shows baseline values of the biological variables, birth weight and weaning weight.
Mean evolution over time of the biological variables is shown in Fig. 1. A global effect of time was observed for 14 out of the 15 biological variables (F D R0.05). However, when considered at each time point individually, only 9 variables
presented a level significantly different from their basal level for at least one
time point (t=+1,+4 or +24).

Fig. 1. Mean evolution of the biological variables overtime. (*): time measurement at which
the expression of the variable is significantly different (F D R0.05) from the expression at t=0. Vertical bars correspond to + and – SEM at each time point

Table 1. Reference values (at t=0) for the biological variables, birth weight and weaning weight (n=120)

As expected 23], ACTH induced a strong cortisol response peaking 1 h after injection (F D R=1.39e^?10
). Plasma cortisol levels at t=+1 were 2.7-fold higher than and highly correlated with basal levels (r²
=0.63,F D R=1.95e^?14
). They were significantly lower than baseline values at t=+4 (F D R2.2e^?16
) and at t=+24 (F D R=1.39e^?10
), as can be expected from the feedback of cortisol on its own secretion. Plasma glucose
levels showed a slight increase at t=+4 after ACTH injection, but FFA levels showed a large increase at t=+1 (×3.21,F D R2.2e^?16
), with a strong variation across animals, and were almost back to basal levels at
t=+4. The values measured at t=+1 were correlated with basal values (r²
=0.42,F D R=2.30e^?06
). Cortisol induces a weak increase in circulating glucose and also potentiates the
effect of other counter-regulatory hormones 24], 25] and increases FFA levels via acute lipolysis 26].

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 27]–29]. 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=+1 and t=+4. These effects are consistent with the literature 30], 31] and result from the redistribution of leucocytes between blood and tissues 32]. Red blood cell related variables (red cells count, hematocrit and hemoglobin levels)
were decreased after injection of ACTH and remained low at t=+24. Platelets were not influenced by injection of ACTH.

Sex did not influence any of the variables (F D R0.05).

Between-subject variability at t=0

Results of the PCA at t=0 on variables responding to ACTH and birth and weaning weights are shown in Fig.
2. As red cells count (RC), hematocrit (Hct) and hemoglobin (Hgb) 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 (i.e., to the first axes of the PCA) in the dataset. The first dimension shows an opposition
between the proportion of granulocytes (positively correlated with this axis) and
the proportion of lymphocytes (negatively correlated with this axis). The second axis
shows an opposition between cortisol (positively correlated with this axis) and birth
and (to a lesser extent) weaning weights (negatively correlated with this axis) 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. 2. PCA on the biological variables identified as responding to ACTH, the birth and weaning
weights at t=0. Colors symbolize the sex: Black = Female; Red = Male a Projection of the individuals on dimensions 1 and 2; b Projection of the variables on dimensions 1 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
(see Fig. 3). The first component of the multi-level PCA opposes the observations at t=0 and t=+24 (positive coordinates on this axis) to the observations at t=+1 (negative coordinates on this axis), this time step corresponding to the peak
of cortisol. The second component opposes the observations at t=+4 (positive coordinates on this axis) to the observations at t=+1 (negative coordinates on this axis). 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 (high measures at t=+1), on one side, and lymphocytes and monocytes (high measures at t=0/+24), on the other side. The second axis shows an opposition between glucose (positively
correlated with this axis, high measures at t=+4) and cortisol, FFA and red count (negatively correlated with this axis, high measures
at t=+1).

Fig. 3. 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=+1 and other variables at t=0 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
2. Only glucose and FFA levels at t=0 were significantly positively correlated with the level of cortisol at t=+1 (F D R0.05). Correlations between cortisol at t=+1 and other variables at t=+1,t=+4 or t=+24 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=+1 were positively correlated with cortisol at t=+1 and glucose at t=+1,t=+4 and t=+24 was negatively correlated with cortisol at t=+1 (F D R0.05). No other variable was found to be significantly linked to the intensity of
the cortisol level in response to ACTH.

Table 2. Correlation coefficients between the biological variables at t=0,t=+1,t=+4 and t=+24 and cortisol at t=+1 (n = 120)

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 4 time points (t?{0,+1,+4,+24}). Differential analysis revealed 158 DE transcripts (adjusted P0.05) matching 65 unique genes (The complete list with features is provided in Additional
file 1). Among them, 23 genes were differentially expressed at t=+1 (5 down regulated/18 up-regulated), 25 were differentially expressed at t=+4 (8 down-regulated/17 up-regulated) and 17 were differentially expressed at t=+24 (all down-regulated). The only gene DE at both t=+1 and t=+4 was SUCNR1 (Table 3). The adjusted P-values were smaller for tests between t=0 and t=+1 and between t=0 and t=+4 than between t=0 and t=+24 (see Additional file 2). This shows that the transcripts were more differentially expressed between t=0 and t=+1 and between t=0 and t=+4 than between t=0 and t=+24. Main effects of cortisol released by ACTH injection on gene expressions are
thus observed at t=+1 and t=+4 with a return to baseline levels at t=+24.

Table 3. List of 65 unique genes differentially expressed in response to ACTH in pigs (n=30)

HAC performed on the within-subject deviation matrix with the list of DE genes identified
4 groups of genes. Figure 4 shows that the 65 unique DE genes allow for an almost perfect classification of the
observations with respect to their time of measurement. For every cluster, Fig. 5 shows the average evolution of each gene (gray) and the average evolution over the
genes in the cluster (red).

Fig. 4. Hierarchical ascending classification of the 65 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=+24. Numbers on the column dendrogram identify each cluster

Fig. 5. Average evolution of the genes in each of the cluster identified by the HAC on the
65 unique DE genes. Evolution of each gene is translated so that it is equal to 0
at t=0; Gray: Average evolution of each of the genes in the cluster. Red: Average evolution
over all genes in the cluster (cluster 1: 18; cluster 2: 17; cluster 3: 8; cluster
4: 22)

Each cluster was then subjected to a functional analysis (results shown in Additional
file 3). In each cluster, genes were DE (F D R0.01) at each time point except for t=+24 in cluster 3 (F D R=0.57).

The first cluster (17 genes) was characterized by genes increasing with a peak of
expression at t=1 and stable between t=+4 and t=+24. 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 (FKBP5, IL7R, CEBPD, CEBPB and NFKB1A in cluster 4). FKBP5 (FK506 binding protein 51) is a decisive factor for the physiological stress response
33] and has an important role in stress-related phenotypes 34]. It modifies steroid hormone receptor sensitivity 35]. CEBPB, DUSP1, FKBP5 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 (C/EBPs) are key regulators of
cell differentiation and are also involved in the expression and production of inflammatory
cytokines 36]. The increase of Period 1 gene (PER1) expression in peripheral blood cells by glucocorticoids was previously reported
in humans 37]. Physical and inflammatory stressors induce the release of the adrenal glucocorticoid
hormone that rapidly alter the expression of PER1 in peripheral tissues through a GRE enhancer present in the gene promotor 38]–40]. This gene is involved in the circadian rythm, in which the glucocorticoid mechanism
plays a predominant role 41]. Another DE gene DDIT4 (regulated in development and DNA damage response 1) was described as a surrogate
biomarker of the efficiency of glucocorticoid receptor blockade in skeletal muscle
42]. Britto and collaborators showed that DDIT4 expression was low under basal conditions but was highly increased in response to
several catabolic stressors, like hypoxia and glucocorticoids 43]. Glucocorticoids were shown to up-regulate DUSP1 in peripheral tissues 44] but constrain the increase of DUSP1 gene expression in the central components of the HPA axis 45]. In vitro studies have shown that glucocorticoid suppression of some MAP-kinase dependent
cellular processes depends on glucocorticoid mediated up-regulation of DUSP1 gene expression 46].

The second cluster (17 genes) was characterized by genes with an increase between
t=0 and t=+4 and a decrease between t=+4 and t=+24. This cluster with genes up-regulated at t=+4 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. CD14 gene is a component of the innate immune system and has been shown to be sensitive
to stress in pigs 47]. MEGF9 gene was shown to be induced by cortisol in human fetal cells in vitro 48].

The third cluster (8 genes) includes the genes decreasing between t=0 and t=+4 and returning to a basal level between t=+4 and t=+24. No ontology was significantly enriched by genes of this cluster. It is interesting
to underline here the ALOX15 gene (arachidonate 15-lipoxygenase) 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 49].

The genes related to the fourth cluster (22 genes) decrease between t=0 and t=+1, increase between t=+1 and t=+4 and decrease between t=+4 and t=+24. The fourth cluster corresponds to genes with an overall expression decreasing
between t=0 and t=+24. They are significantly linked to biological processes such as protein phosphorylation
and kinase activity. Among the genes involved in this cluster ARHGAP31 and ARHGAP family genes were found to be differentially expressed in macrophages treated with
dexamethasone 50], 51].

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” (F D R0.05). These genes are thus differentially expressed over time even when adjusting
for changes in white cells populations. Among them, 34 genes had their expression
significantly negatively influenced by L/G ratio (see Additional file 1), 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 2, over-expressed
at t=+4 (17/17) and cluster 1, over-expressed at t=+1 (12/18) and to a lesser extent as genes of cluster 4, under-expressed at t=+24 (5/22). No gene of cluster 3 was significantly explained by L/G ratio. Results
of the functional analysis of this list of 34 genes are shown in Additional file 4. 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 65 DE genes and the biological variables were then subjected to a multilevel PLS.
Figure 6 shows that the first axis of the multilevel PLS opposes observations at t=+1 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. 6. 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=+24. b projection of the variables; Blue: biological variables; Red: gene expressions; 10
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, RGS2, RHOB, PER1, FKBP5, CEBPD, DDIT4, CPT1A and DUSP1. 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 SUCNR1, SLCO2B1, FBP1 and LOC396700. These genes belong to the third cluster and are related to glycolysis and glycogenesis.
SUCNR1 (succinate receptor 1) is decreased at t=+1 and increased at t=+4. Succinate has a wide range of metabolic actions and regulates the functions of
macrophages 52]. The axis is negatively correlated with CD14, CLC4D, CHIT1, MEGF9 and C2H19orf59. 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 (DDIT4, DUSP1, FKBP5, IL7R, NFKBIA, PER1, RGS2, RHOB, Fig. 7) are functionally connected to each other by NR3C1. The NR3C1 (nuclear receptor subfamily
3, group C, member 1) 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 9]–13]. Mutations in NR3C1 have been previously demonstrated to be associated with generalized glucocorticoid
resistance 53]. 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. 7. Gene network related to glucocorticoid response in whole blood transcriptome 1 h after
ACTH injection This network corresponds to the genes up-regulated 1 h after ACTH injection
(cluster 1, green nodes). It combines bibliographic (best enrichment score network
=45) and regulatory relationships (genes co-regulated by the same regulator in blue
with the highest enrichment score (p-value=2.00E?09), green lines) proposed by Ingenuity software. Cortisol has the highest plasma
level at t=+1 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 1 and shown
as ACTH responsive were found differentially expressed in stress-related investigations.
Five genes found in our study (CXCR4, DUSP1, FKBP5, IL7R, TXNIP) were proposed as markers of differential glucocorticoid sensitivity 54], 55]. NFKBIA, DUSP1, CEBPD, FKBP5 genes were also found to be associated with up- and down-regulated clusters in response
to continuous 24 h cortisol infusion 56]. Ponsuksili and collaborators 57] describe NFKBIA, CEBPB and CEBPD as genes of which hepatic expression levels are correlated with plasma cortisol concentrations.
Up-regulation of PER1 gene upon GR activation was confirmed by genome-wide study of glucocorticoid receptor
binding sites in neuronal PC12 cells 58]. However, DDIT4 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=0,t=+1,t=+4 or t=+24 were significantly correlated with the level of cortisol at t=+1, only two genes were identified: TRMT2A (F D R=0.04), a gene involved in the methylation of tRNA, and LOC100626276 (F D R=0.04), a gene of which function has not been identified yet. There is a negative
relationship between the expression of cortisol at t=+1 and the expression of these two genes at t=0 (?0.45 and ?0.63 for TRMT2A and LOC100626276, respectively). This implies that when their baseline expression is higher, the intensity
of the cortisol response to ACTH decreases.