Transcriptome profiles of hypothalamus and adrenal gland linked to haplotype related to coping behavior in pigs

The hypothalamic-pituitary-adrenal (HPA) axis is an important component of neuroendocrine stress regulation and coping behavior. Transcriptome profiles of the hypothalamus and adrenal gland were assessed to identify molecular pathways and candidate genes for coping behavior in pigs. Ten each of high- (HR) and low- (LR) reactive pigs (n = 20) were selected for expression profiling based haplotype information of a prominent QTL-region on SSC12 discovered in our previous genome-wide association study (GWAS) on coping behavior. Comparing the HR and LR pigs showed 692 differentially expressed genes (DEGs) in the adrenal gland and 853 DEGs in the hypothalamus, respectively. Interestingly, 47% (17 out of 36) of DEGs found in both tissues were located in GWAS regions identified on SSC12, indicating that there are significant functional positional candidate genes for coping behaviour. Pathway analysis assigned DEGs to glucocorticoid receptor signaling in the adrenal gland. Furthermore, oxidative phosphorylation, mitochondrial dysfunction, and NGF signaling as well as cholecystokinin/Gastrin-mediated were identified in the hypothalamus. We narrowed the list of candidate genes in GWAS regions by analyzing their DEGs in the HPA axis. The top identified transcripts, including ATP1B2, AURKB, MPDU1 and NDEL1 provide evidence for molecular correlates of coping behavior in GWAS regions.


Results
Pigs were selected for transcriptome profiling of hypothalamus and adrenal gland tissue based on the backtest haplotype and phenotype. We first calculated the haplotype association, testing the markers ALGA0066975, ALGA0121951, ASGA0055092, ASGA0105202, H3GA0034753, and MARC0073387 in all 294 individuals. All candidate markers were located on SSC12 position 55-56 Mb (Ensembl_Sscrofa_10.2) and were previously reported as single markers associated with backtest traits 22 . In total, 4 haplotype variants of 294 piglets were used to test the association against the backtest phenotype. Haplotype 1 (A/G/A/C/A/A, n = 180) and 2 (G/A/G/A/ G/C, n = 14) were found to be associated with backtest traits (Table 1). Based on this information together with backtest parameters, animals were divided into two groups, the high reactive (HR) group and low reactive (LR) group, each including 10 animals. A summary of phenotypic data such as the backtest characteristics of latency, duration and frequency, and cortisol levels is shown in Table 2. HR and LR pigs clearly showed significant differences in the backtest parameters for latency (p = 0.0002), duration (p < 0.0001) and frequency (p < 0.0001). Latency, a measure of time elapsed to first struggle, was nearly five times higher in the LR than the HR pigs. In contrast, cumulative struggling duration of the HR pigs was about four times higher than that of the LR group. The HR pigs showed a higher frequency of struggling attempts than the LR pigs during the 1-minute testing period. Cortisol levels did not differ significantly between the LR and HR pigs.
Differentially expressed genes and pathway analysis. Hypothalamus. From a total of 47 880 probe sets on the snowball array, 17 351 passed filtering steps and were ready for further analysis. Mixed model analysis  Table 1. Statistical testing of association between 4 haplotype variants (H1-H4) of 294 piglets against backtest parameters duration, frequency and latency. Data set lists F-statistics and associated probabilities (Prob > F) for each of the estimated haplotypes. Rare haplotypes are defined as haplotypes with a frequency lower than 0.05.
Adrenal gland. In the adrenal gland, 19887 out of 47880 probe sets remained after filtering for downstream analysis. A mixed-model analysis revealed 692 DEGs between the two groups at p < 0.05 (Supplementary Table S2). Of these, 512 DEGs were annotated by IPA. Three hundred sixty-three DEGs were down-regulated in the HR compared to LR pigs. Twenty-three DEGs exhibited an absolute fold change greater than 1.5. The top five DEGs with increased mRNA transcript abundance in HR compared to LR pigs were IGSF5, MLKL, EDRF1, MGAM, and PCDH7. The top five genes with lower expression levels in the HR pigs were ATP1B2, OAS1, PTGFRN, PRLR and SEPW1 (Table 5). The ATPase Na + /K + transporting subunit beta 2 (encoded by ATP1B2) showed the highest absolute fold change (FC = −4.04, p < 0.0001) with expression levels in HR pigs lower than levels in LR pigs. Up-regulated (HR > LR) genes were mainly associated with PPARα/RXRα activation, whereas down-regulated (HR < LR) genes were related to EIF2 signaling, and glucocorticoid receptor signaling, estrogen receptor signaling (Table 6). Tissue-wide analysis of differentially expressed genes of both tissues in GWAS regions: In order to identify candidate genes in GWAS regions of both tissues, mixed-model analysis with yielded 1428 differentially expressed genes between the LR and the HR group in the tissue-wide analysis of both hypothalamus and adrenal gland. To functionally characterize the genetics of coping behavior and identify candidate genes, we overlaid these identified DEGs with our previous GWAS results 22 , as shown in Table 7. A total of 36 out of 1428 DEGs were found in proximity (50 kb upstream and downstream) of trait-associated SNPs. The major share (47.22%) was

Discussion
Improving animal welfare and reducing the psychosocial stress of domestic pigs is a vital issue in modern pig production. Several studies explored coping behavior from various angles, focusing on behavioral 3 , physiological 23 and immunological 24 aspects. Krause et. al demonstrated coping style-dependent central physiological differences in autonomic nervous system regulation, which is controlled by the hypothalamus 25 . For our study, we selected animals based on the haplotype of candidate markers previously found to be highly associated with coping behavior 22 . These genomic regions, harboring several genes, exhibited high linkage disequilibrium (LD) 22 . The above mentioned haplotypes may be in close LD with the causal variant or may themselves consist of the causal variants of interest 26 . In addition, the biological functions of transcripts of specific tissues involved in the hypothalamic-pituitary-adrenal (HPA) axis from the same animals were taken into account. Coping behavior associated haplotype information in synergy with the evaluation of biological transcript function facilitates the discovery and broadens the knowledge of causal variants for complex phenotypes like coping behavior and thus might provide promising new biomarkers for genotype-based selection and animal welfare. Significant differences between the HR and LR groups in all three backtest parameters (latency, duration, and frequency) emphasize an association of genetic based early-life coping style with differential mRNA expression in our study. Several previous transcriptomic studies focused on the differential mRNA expression associated with adrenal sensitivity to ACTH and psychosocial stress [27][28][29][30] . The present study is the first to report associations between genetics, temperament, and differential gene expression in porcine hypothalamus and adrenal gland tissues. Preslaughter handling including the use of electronarcosis to sacrifice the animals may affect animal welfare and stress response, particularly in brain regions. Animals of both groups were treated the same way. Plasma glucocorticoid levels are considered good markers of stress 31 , but in this study, the total plasma cortisol level measured by ELISA showed no significant difference between the HR and LR groups; although, it tended to be higher in the LR group. Our previous study measured cortisol levels in the absence of any backtest experiments in 475  www.nature.com/scientificreports www.nature.com/scientificreports/ individual pigs and showed levels of 93.9 ± 34.6 ng/ml (mean ± SD) 32 . The range of this baseline measurement lies between the cortisol measurement ranges of the HR and the LR group.
This observation suggested that the HR and LR pigs were likely not exposed to stressful conditions prior to sample collection and that the gene expression pattern was not directly influenced by the cortisol level. In this study the samples were taken at adult stage. A prior study has shown that repeated backtests reflect personality and coping strategy with a moderate intra-individual consistency and heritability 19 . More than 3000 animals, of which the animals used here represent a subset, it was previously shown that the correlations among the backtest at different time points (ages) (rs = |0.19-0.44|) strongly supports the backtest traits as reliable parameters of personality and coping style and for long life and inheritance 19 .
Differentially-expressed transcripts in the hypothalamus. Due to its pivotal role in stress response, the transcriptomic signature of the hypothalamus is of particular interest and may shed light on the molecular basis of coping-related behavioral traits. Altered mRNA transcript abundances in the HR compared to LR group in the hypothalamus were associated with signaling of cholecystokinin (CCK), nerve growth factor (NGF), and their interaction as well as with oxidative phosphorylation, oxidative stress and mitochondrial dysfunction.
In the hypothalamus, we found that Cholecystokinin/Gastrin-mediated signaling is upregulated in the HR pigs. Besides its well-known role as a peptide hormone within the gastrointestinal system, cholecystokinin is found extensively throughout the central nervous system (CNS) 33,34 . Previous studies have reported a wide variety of behavioral and autonomic actions in different species, such as the suppression of exploratory behavior 35 , a modulatory role in memory formation and learning 36 , anti-opioid activity [37][38][39] , anxiety induction 40 , thermoregulation, and involvement with the dopaminergic reward system 41 . Human as well as rodent studies clearly suggest a positive causal relationship between CCK and anxiety 42 . A robust and dose-dependent release of ACTH and cortisol could be observed after intravenous injection of the CCK-B receptor agonist pentagastrin 43,44 . Preliminary data showed a relative resistance to cortisol feedback inhibition 43 . Consequently, CCK could be involved in an acute activation of the HPA axis, which is independent from the actual plasma cortisol level and in this way promotes its activity even under conditions of already elevated cortisol. It is noteworthy to mention that the anxiogenic and HPA effects of CCK are independent phenomena being coordinated by multiple pathways [44][45][46] . Beyond that and strikingly, CCK-8 has been implicated in inhibiting oxidative stress-induced neurotoxicity via an anti-oxidative stress pathway 47 .
Higher basal gene expression of transcripts associated with the IPA biofunction NGF signaling in the HR group were observed in our study. As a highly conserved neurotrophin, nerve growth factor plays a pivotal role in the survival and maintenance of sympathetic and sensory neurons 48 . Interestingly, a relationship between CCK and NGF has been reported previously: intraperitoneally injected CCK-8 was shown to stimulate NGF synthesis via activation of CCK receptors 49 . Taking the role of CCK-8 regarding oxidative stress response into account, activation of CCK-/gastrin mediated signaling can be seen as one step in an oxidative stress-induced neuronal recovery process carried out by NGF. Animal and human studies report reduced body weight, increased anxiety,  www.nature.com/scientificreports www.nature.com/scientificreports/ and mood changes after NGF administration 50,51 . Furthermore, NGF injection into the brains of cholinergic function-compromised rats induced hyperactivity and fear 49 .
The gene coding for the highly conserved transcription factor SRF presented with lower expression levels in the HR compared to LR pigs. Many immediate early genes such as c-fos are regulated by SRF. A regulatory function of SRF signaling is mediated by MICAL2, an atypical actin-regulatory protein 52 . In our study, MICAL2 expression was elevated in the hypothalami but decreased in the adrenal glands of HR pigs compared to LR pigs. A basic capability of MICALs is the generation of redox potential, via their mono-oxygenase domains, and in turn reactive oxygen species (ROS) 53 . For this reason, we speculate that hypothalami of LR individuals face heightened ROS production due to increased MICAL2 expression. Earlier studies have emphasized the role of MICAL2 and several adjacent genes in mediating anxiety 54 .   www.nature.com/scientificreports www.nature.com/scientificreports/ MUC1 was one of the top ten hypothalamic DEGs in our study (FC = −1.54, p < 0.0014). Its cytoplasmic tail is able to bind TP53, increasingly under conditions of genotoxic stress 55 . TP53, also one of the hypothalamic DEGs, aligns with QTL regions on SSC12 (55.2 Mb) and showed decreased transcript abundance in the HR pigs. Both molecules, MUC1 und TP53 are associated with the TP53 response element p21 gene promotor, leading to an activation of p21 and finally cell cycle arrest 55 . Additionally, MUC1 expression increases as a cellular response to oxidative stress 56 ; hence, its elevated expression in LR compared to HR pigs can be interpreted as an indicator www.nature.com/scientificreports www.nature.com/scientificreports/ of a derailed reactive oxygen species (ROS) balance, which may be caused by increased MICAL2 expression in the LR pigs.
Due to their crucial role in cellular physiology, mitochondria are particularly sensitive to changes in cell homeostasis. It is therefore not surprising that we found different activation patterns of pathways related to oxidative phosphorylation and mitochondrial dysfunction in the LR and HR groups. DEGs associated with mitochondrial complexes I and III were lower in HR individuals. Even partial inhibition of complex I can cause enhanced ROS production, leading over time to oxidative stress and eventually neuronal damage 57 . Many psychiatric and www.nature.com/scientificreports www.nature.com/scientificreports/ neurological disorders are correlated with compromised mitochondrial function 58 , and increased anxiety-like behaviors have been observed in mice with mitochondrial dysfunction 59 . In a zebrafish study investigating chronic unpredictable stress (CSU)-induced anxiety, brain proteome profiling revealed differentially regulated proteins predominantly involved in mitochondrial function, glycolysis, oxidative stress, and the hypoxia stress response, emphasizing that mitochondrial dysfunction may be causal for the anxious phenotype 60 .
Differentially-expressed transcripts in adrenal gland. As part of the HPA axis, the adrenal gland produces glucocorticoids which are able to affect almost any vertebrate cell by binding to the glucocorticoid receptor. Moreover the gland secretes catecholamines into the circulatory system in the course of sympathetic activation. These physiological roles make this tissue a promising target for elucidating the molecular pathways of coping behavior. The DEGs in adrenal gland have higher q values, despite using samples from the same animals we used for hypothalamic differential expression analysis. The adrenal gland constitutes the last component of the HPA axis which might explain the limited impact on differential gene expression between the HR and the LR pigs. Nevertheless, this organ is highly responsive to exogenous stressors and in fact, we found that the most significant adrenal signatures differing between the LR and the HR pigs were associated with glucocorticoid receptor signaling.
Cortisol and other glucocorticoids bind to the glucocorticoid receptor to execute gene regulatory functions in development, metabolism, and immune response. IPA associated 10 DEGs in our present study to the glucocorticoid receptor signaling pathway, including BAG1, BCL2L1, CD247, ELK1, NFKBIA, POLR2G, POLR2H, SMAD2, STAT5B, TAF3 and HDAC6. Overall, activation of this pathway was reduced in HR compared to LR pigs. Downregulation of glucocorticoid receptor signaling ultimately entails downregulation of NFKBIA in the HR group. Expression of the NFKBIA gene, encoding for an inhibitor of NFKB-mediated inflammation, is stimulated by glucocorticoids 61,62 and was shown to be strongly correlated with plasma cortisol concentrations 32 .

Differentially expressed genes of both tissues in GWAS regions.
Approximately 47% (17 out of 36 genes) of the DEGs in both tissues were located in GWAS regions which have been identified on SSC12 22 . These candidate genes in combination with their biological function will be selected for further analysis of the causal relationship of SNPs and behavioral traits. Results from the previous GWAS linking backtest parameters and coping style with significantly associated markers revealed a cluster of SNPs with LD in the region between 54 and 56 Mb of porcine chromosome SSC12 using the pig genome assembly Sscrofa10.2 22 . The corresponding region in the new pig genome assembly (Sscrofa11.1) is located between 52 and 54 Mb on SSC12. Among the tissue-wide DEGs in GWAS regions, three highly significant genes (FDR < 0.0001) located in the GWAS region on SSC12 deserve particular attention: ATP1B2, MPDU1 and AURKB. ATP1B2 exhibited significantly lower expression (FDR < 0.0001) in the LR group. By using qPCR for validation, we confirmed both differential expression of ATP1B2 in the adrenal gland (p < 0.0001, FC = −4.04) and in the hypothalamus (p = 0.0175, FC = −1.54). Consistent with this, a study with two substrains of Wistar Kyoto rats (high mobility and low mobility), which were subjected to forced swim test, revealed differential expression of ATP1B2 in their hippocampi 63 . ATP1B2 (also known as AMOG) encodes the beta 2 subunit of Na + /K + -ATPase. This sodium pump is responsible for the electrochemical gradient of Na + and K + ions and thus maintains membrane potential and cellular homeostasis 64 . In addition, dysfunction of AMOG might influence the ionic and osmotic regulation, contributing to neuronal hyperexcitability 65 . ATP1B2 −/− mice showed rapidly worsening motor impairment and degeneration of neural cells in the central nervous system 66 . A recent study suggest that the neurohistopathological changes seen in spongy degeneration with cerebral ataxia (SDCA) in Malinois dogs might be caused by a SINE insertion into ATP1B2, leading to aberrant RNA splicing and reduced ATP1B2 protein expression 67 . Mannose-P-dolichol utilization defect 1 protein, encoded by the MPDU1 gene, plays an integral role in the synthesis of glycosylphosphatidylinositols and lipid-linked oligosaccharides 68 . Homozygous point mutations in the MPDU1 gene have been associated with congenital disorder of glycosylation type If in humans, a disease featured by, inter alia, psychomotor retardation and seizures 69 . The underlying cause of seizures is seen in abnormalities of neuronal activity in the brain due to excessive or synchronous firing 70 . The link between MPDU1 and the occurrence of seizures in humans makes it tempting to speculate, that this gene participates in one way or another in neuronal activity regulation and that lower MPDU1 expression could lead to abnormally excessive neuronal firing. AURKB encodes a member of the aurora kinase subfamily of serine/threonine kinases. Recently, it has been suggested that AurkB is involved in the regulation of neuronal development and axonal outgrowth 71 . Experiments, using pharmacological and genetic methods to impair AurkB function rendered motor axon morphology to a truncated and abnormal phenotype, whereas elongated axonal outgrowth has been observed after AURKB overexpression 71 . Both AURKB and MPDU1 were found to be highly differentially expressed between HR and LR of both tissues but the functions of these genes with regard to coping behaviour still remain unknown. Another gene in this region, NDEL1, plays a role in multiple processes including cell signalling, neuronal migration and neurite outgrowth. Ndel1 is responsible for motor protein (dynein) regulation together with LIS1. Disruption of this complex leads to impaired neuronal migration by altering nuclear-centrosome coupling 72 . By its ability to bind DISC1, the NDEL1/Lis1 complex is involved in the development of the neurobehavioral symptoms of schizophrenia 73 . A genome-wide study investigating the genetic factors of schizophrenia found an association between the enzyme activity of NDEL1 and disease symptoms, with lower levels of NDEL1 activity in the case group 74 . A well described human neuronal migration disorder caused by aberrant Lis1/Ndel1 function is Miller-Dieker syndrome (MDS) 75 . Affected individuals suffer from severe mental retardation as well as epileptic seizures, suggesting that NDEL1, like MPDU1, might contribute to excessive neuronal firing. Located on SSC7, FAH belongs to the promising DEGs in QTL regions of both tissues. Fumarylacetoacetate hydrolase (FAH) catalyzes the last step of tyrosine degradation 76 , and tyrosine itself is the amino acid precursor of catecholamine synthesis 77  www.nature.com/scientificreports www.nature.com/scientificreports/ elevated mRNA transcript abundance of FAH in the LR group observed in our study suggests increased tyrosine degradation and thus an abated supply of catecholamine precursors with corresponding effects on the sympathetic activity of this group. In contrast, absence of FAH expression in mice brains produced a phenotype with markedly increased myelination and altered social behavior 78 .
In summary, transcriptomic profiling in combination with information on genetic variance and backtest behavior proved to be a versatile approach for understanding genotype-phenotype-mapping associated with targets of hypothalamic and adrenal gland gene expression and the physiological processes related to backtest behavioral traits. Our study demonstrated haplotype associated, coping pattern linked transcripts in the hypothalamus and the adrenal gland. These findings revealed meaningful biological pathways like oxidative phosphorylation, mitochondrial dysfunction, cholecystokinin/gastrin-mediated signaling, and NGF signaling -the latter two working hand in hand to counteract oxidative stress-induced neurodegenerative effects in the hypothalamus and themselves affecting behavioral traits. The glucocorticoid receptor pathway is plausible molecular pathway that plays a role in this context in the adrenal gland. 47% of the DEGs of both tissues were located in GWAS regions which have been identified on SSC12. Among them, ATP1B2, MPDU1, AURKB and NDEL1 represent the most promising candidates for further analyses. The gene ATP1B2 is involved in the maintenance of electrochemical gradients and neuronal excitability, NDEL1 is associated with both neuronal migration and neurite outgrowth and despite clear evidence for a functional role in the nervous system; the exact link of MPDU1 and AURKB to behavioral traits still needs to be unraveled. Altogether, the results of our study increase our understanding of causal variants for complex phenotypes like coping behavior in pigs. However, it is advisable to confirm the role of the identified genes and pathways for coping behavior in other pig populations and other ages of pigs, as alternative pathways may be relevant in other populations and ages. By narrowing down the list of candidate genes in GWAS regions, this study provides evidence for molecular correlates of coping behavior, which may have the welfare-enhancing potential to be used as biomarkers for animal wellbeing or genotype-based selection.

Material and Methods
Ethics statement. Animal care and tissue collection procedures followed the guidelines of the German Coping behavior and haplotype estimation. Animals were provided by the Leibniz Institute for Farm Animal Biology (FBN, Dummerstorf, Germany). 294 German Landrace pigs used in this study were subjected to backtests on days 5, 12, 19 and 25 post natum which was previously reported 22 . Briefly, each individual piglet was swiftly turned on its back and put onto a V-shaped device padded with cellulose. To avoid escape, the piglets were held manually in that position with slight pressure, adapted to the intensity of their movements. The test commenced immediately after reaching immobility and lasted for 60 s. Latency (L), the time elapsed until the first struggling attempt, duration, the cumulative struggling time (D) within the 60 s test, and frequency (F), the number of struggling attempts, were recorded. Moreover, total latency (tL), total duration (tD), and total frequency (tF) were calculated for each animal by summing up the relevant parameters at all ages tested.
Based on our previously reported of a prominent QTL region harboring a coping-behavior associated haplotype block on SSC12, we used the following 6 markers (ALGA0066975, ALGA0121951, ASGA0055092, ASGA0105202, H3GA0034753 and MARC0073387), located in a linkage disequilibrium (LD) block within the GWAS region of SSC12 ( Supplementary Fig. S1), for haplotype estimation and trend regression analysis 22 . To estimate the unobserved haplotype frequencies, the haplotype estimation process in JMP Genomics (SAS Institute, Cary, NC, USA) was used. This process invokes an expectation-maximization (EM) algorithm to estimate haplotype frequencies. Estimates of haplotype frequencies were used as input for the Haplotype Trend Regression process, in order to further determine the particular haplotype from the 6 SNPs in LD regions. To test for association of each haplotype with backtest parameters, PROC GLIMMIX was used with sex as a fixed effect and sire as a random effect. In total, 4 haplotype variants of 294 piglets were used to test the association against backtest phenotype. The data set lists the F-statistics and associated probabilities for each of the estimated haplotypes. HR animals are homozygous carriers of haplotype 1 (A_G_A_C_A_A), whereas all LR animals are carriers of haplotype 2 (G_A_G_A_G_C). A total of 20 pigs with high-reactive (HR, n = 10) or low-reactive (LR, n = 10) coping patterns together with the information of the haplotype were selected from a pool of 294 German Landrace (DL) piglets.
Sample collection. Pigs were weighed and slaughtered by electronarcosis followed by exsanguination. The average age of the 20 animals was 157 days post natum. Blood samples were collected in 50 ml Falcon tubes containing 1 ml of 0.5 M EDTA and immediately stored on ice. Subsequently, plasma sample separation was performed, and samples were kept at a temperature of −80 °C until use. An enzyme-linked immunosorbent assay (ELISA, DRG, Marburg, Germany) was carried out in duplicate as described previously 22 in order to provide total plasma cortisol levels. Additionally, tissue samples for the determination of transcript abundances were taken immediately after exsanguination. For the dissection of hypothalamus, a stereotaxic atlas of the porcine brain was used as reference guide 79 . The whole hypothalamic area, including the paraventricular nucleus, was extracted surgically out of the left and right hemispheres after quickly harvesting the brains. In brief, the skull was opened www.nature.com/scientificreports www.nature.com/scientificreports/ with a saw. The optic chiasma and the nucleus mamillaris anterior (Mammilary body) up to the thalamus were removed. Laterally of the anterior commissure along the thalamus was cut. The collected tissue corresponded to the width of the mammilary body. This procedure has ensured that the Nucleus paraventricularis hypothalami and the Nucleus arcuatus hypothalami were part of the collected tissue. Middle portions of both adrenal glands were excised. After taking tissue from the hypothalamus and adrenal gland, the samples were snap frozen in liquid nitrogen and stored at −80 °C until use.

RNA isolation target preparation and hybridization.
Total RNA was isolated from hypothalamus and adrenal gland tissue samples using the TRI Reagent (Sigma-Aldrich, Taufkirchen, Germany) and RNeasy kit (Qiagen, Hilden, Germany). The RNA was quantified using the NanoDrop ND-1000 spectrophotometer (Peqlab, Erlangen, Germany), checked for integrity by performing agarose gel electrophoresis (1% agarose gel), and stored at −80 °C until use. To prepare the samples for microarray analysis, amplified sense-strand cDNA was generated using the Ambion WT Expression Kit (Ambion, Austin, TX, USA). Next, the cDNA was fragmented and biotin-labeled using the Affymetrix GeneChip WT Terminal Labeling Kit (Affymetrix, Santa Clara, CA, USA). Each individual sample was hybridized on a genome-wide snowball array (Affymetrix, Santa Clara, CA, USA), containing 47880 probe-sets. After a staining and washing step, the arrays were scanned and processed using Affymetrix GCOS 1.1.1 software.
Microarray data processing. The data was pre-processed using Affymetrix Expression Console 1.3.1.187 software (Affymetrix), and the normalization was done using the RMA (robust multichip average) algorithm. The DABG (detection above background) algorithm was used to filter present (expressed) genes. Probes-sets present in less than 80% of the total number of samples and those representing miRNAs were excluded from further analysis. Differential gene expression analysis between the LR and HR group was assessed by running mixed model procedure in JMP Genomics 7.0 (SAS Institute, Cary, NC, USA) with sex and coping style as fixed effects and sire as a random effect. The adjusting for multiple comparisons across the Type 3 tests for all of the fixed effects was calculated using the post hoc Tukey-Kramer test. To control for multiple testing, p-values were converted to a set of q-values 80 . The level of significance was set at p < 0.05 which corresponds to q < 0.17 for hypothalamus and q < 0.47 for adrenal gland, respectively. Differentially expressed genes (DEGs) were used for pathway analysis using Ingenuity Pathway Analysis (IPA) (Ingenuity Systems, Redwood City, CA). IPA categorizes genes based on annotated gene functions and statistically tests for over-representation of functional terms within the gene list using Fisher's exact test.
The identification of candidate genes with differential expression between the LR and the HR group from both tissues was implemented using the mixed model procedure in JMP Genomics 7.0 (SAS Institute, Cary, NC, USA) with sex, coping style and tissue as fixed effects and sire as a random effect. The model was combined with a repeated statement with animal as the subject effect, in order to take into account the correlations among measurements made on the same animal in both tissues. The post hoc Tukey-Kramer method was used for multiple comparison adjustments. To control for multiple testing, the FDR was set to 0.05 81 .
Real-time quantitative PCR (qPCR) validation. DEGs were selected for qPCR validation based on their functional roles in the IPA significant pathways using a 48 × 48 Dynamic array with an integrated fluidic circuit (IFC) on the BioMark HD Real-time PCR System (Fluidigm, South San Francisco, CA, USA). cDNA was generated using 2 µg of total RNA, Superscript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA), and oligo(dT), with specific target amplification (STA) and Exonuclease I treatment. A total of 2.25 µL of the sample treated with STA and Exo-I, 2.5 µL of SsoFast EvaGreen supermix with low ROX (Biorad, Hercules, CA, USA), and 0.25 µL DNA binding dye sample loading reagent were loaded for each sample inlet. Assay inlets were loaded with 2.25 µL DNA suspension buffer, 2.5 µL assay loading reagent, and 0.25 µL of a 100 µM primer solution (forward and reverse). Temperature scheme for the qPCR reaction was as follows: An initial denaturation step for 60 s at 95 °C, followed by 30 cycles consisting of denaturation at 95 °C with 5 s each and annealing at 60 °C for 20 s. Reference genes (HPRT1 and RPS11) were selected based on microarray data as invariant genes between LR and HR group. All measurements were performed in duplicate. Primer sequences can be accessed in Supplementary Table S3.

Data Availability
The data used in this study was deposited in the database of the National Center for Biotechnology Information Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo) [GEO: GSE109155, GSM293170-2933189 and GSM2933190-2933209].