Phenotypic assessment of temperament over the years
Young calves, such as at weaning stage, usually show a fearful response to handling due to the novelty of being handled coupled with greater exposure to a variety of potential stressors including breaking the dam-calf bond, change of diet, and re-grouping with unfamiliar animals, as described by Enríquez et al. . Furthermore, studies have shown that the average temperament score decreases as the frequency of handling or additional animal–human interaction increases, as long as the stockperson uses a calm approach to handling the animals [3, 5, 45, 46]. In other words, one would expect the average temperament to decrease (animals become more docile) as the age increases. On the one hand, experimentally-designed datasets are often associated with frequent human–animal interactions. On the other hand, in a seedstock/commercial setting, the handling frequency (i.e., human–animal interactions) decreases as the animals get older. With that, temperament scores taken on mature animals in this dataset were collected at or near the time calves were weaned which could lead to elevated stress responses while handling. Consequently, in practical scenarios an opposite behavioral pattern could be observed: the average temperament scores slightly increased as the animals aged (Fig. 2a), suggesting a sensitization across the years instead of habituation. In short, habituation can be defined as decreased responsiveness to a stimulus while sensitization is an increased responsiveness to a stimulus . Habituation is often associated with adaptation .
We compared younger animals (YT, < 2 years-old) versus older than 2 years (CT, 3 years-old and older), in which differences can be attributed to factors related to the age and the recording event. First, 320 days is the minimum age when animals were recorded for temperament (YT), allowing the calves to acclimate to the challenging post-weaning environment . Second, temperament at ages older than 2 years was recorded on cows around the weaning event. Cows differ in maternal behavior and in their response to weaning of their calf [47, 48], and separation from the calf is likely to trigger different motivations and physiological and behavioral responses during recording of CT compared to YT, although cows had typically been weaned a few weeks before recording of CT. Cows could also change their behavior to a flightier behavior (e.g., scores 2 or 3; Table 2) as a consequence of mothering instincts, in order to be more protective of their progenies.
The average cow at weaning temperament score (animals older than 2 years) increased slightly with age (Fig. 2a). The magnitude of changes of temperament scores between adjacent age groups was slightly greater between YT to CT (increased by 12%) than within CT (maximum increase between adjacent years equal to 2%), which also supports the argument that different lifetime events might evoke different motivational systems and capture context-specific behavioral responses. Furthermore, supporting these results, the phenotypic correlations between YT and CT were weakly positive [0.37; Tables 6 and (see Additional file 1: Table S8)], while strong phenotypic correlations were observed between CT scored in adjacent years and weaker correlations were estimated between larger age-group gaps (from 1.00 to 0.18; Fig. 4b).
Previous negative experiences can trigger subsequent fearful responses to handling due to memory acquisition [9, 16, 17, 49]. Although the number of human–animal interactions is expected to cumulatively increase across the animals’ lifetime, the frequency of interactions per year can decrease (e.g., until yearling the calves are intensively and routinely managed while older animals could be handled less than three times a year). An alternative hypothesis is that animals got, on average, flightier over time due to memorization of cumulative experiences. Such increased reactions to handling can impact the animal and handler welfare as well as impact the behavior of young replacement animals in the herd, because more temperamental dams also raise, on average, more temperamental progeny (Fig. 2b).
Another potential justification for younger animals to be less flighty than older animals is selection. Direct genetic selection based on the official genetic evaluation for temperament would result in a more docile younger population. Based on this theory, flightier animals would be censured on later-in-life temperament measurements. However, in practice, selection is not based on a single trait but rather on a selection index (group of traits) and thus, only animals with extremes breeding values for temperament might have been culled.
Genetic and genomic assessment of temperament over the years
A strong positive genetic correlation was observed between YT and CT (0.84, Table 6) and between pairs of year-groups for CT (average equal to 0.98, Fig. 4). Although these correlations are strong, we observed a similar pattern of phenotypic correlations: stronger within age-groups of CT than between YT and CT. Furthermore, one out of the 24 genes identified for YT by Alvarenga et al.  overlapped with genes located within the top 20 genomic regions based on absolute effect for each age-group for CT (Fig. 8).
The overlapping U6 paralog gene is a noncoding RNA. Alvarenga et al.  systematically reviewed genes that control behavioral indicators in mammalian livestock, and U6 paralogous genes had been previously associated with three behavioral indicators in cattle and two indicators in pigs. Furthermore, a haplotype-based GWAS on yearling temperament in Angus cattle also identified U6 as a candidate gene controlling yearling temperament .
The main overlapping gene between YT and CT, SPTLC2, is located on BTA10 and is associated with neurological diseases in humans  and re-learning mechanisms in rats . The SPTLC2 gene plays a role in the de novo synthesis of sphingolipids (i.e., responsible for signal transduction) by condensing L-serine and palmitoyl-coenzyme A into 3-keto sphinganine, which can be further converted into sphingoid bases . Acid sphingomyelinase (ASM) is a key enzyme in sphingolipid metabolism, and decreased activity of ASM in the dorsal hippocampus has been associated with efficient re-learning in rats . Interestingly, Huston et al.  also reported changes in the ASM enzymes as the rats aged.
Many genes identified for the top 20 SNP-windows based on the absolute value were common across age groups (Fig. 7) and (see Additional file 1: Table S3), which supports the strong genetic correlations observed between CT age groups. Although the genes controlling CT are similar, their effect can be time-dependent. Oliveira et al.  observed an increase and decrease in the magnitude of the SNP effects across time for milk yield and somatic cell score in dairy cattle, respectively. In our study, the magnitude of SNP-window effect increased with age for the top 20 genomic regions (Fig. 6). Furthermore, studies on mental disorders in humans have shown that genes become activated in response to life circumstances, consequently, initiating adverse behaviors . In other words, we hypothesize that as the animal ages, behavioral-related genes are differentially regulated and/or their expression is altered. However, no conclusion can be drawn from the present study. Gene expression studies should be conducted on the genes identified to investigate this hypothesis.
In the same context, the genes that did not completely overlap with all age groups were identified later in life, and in general, they have been previously associated with age-related neural and mental disorders in humans and mice, such as memory and Alzheimer’s disease in older mice. The carboxylesterase 5A gene (CES5A) is located on BTA18, and it was identified as a candidate gene based on the top 20 SNP-regions for cows older than 7 years (see Additional file 1: Table S3). CES5A has been validated as a regulator of male mice fertility [53, 54]. As an expected role of carboxylesterase family genes, the CES5A protein can result in high levels of cholesterol and choline esterase , alterations of which impact cognitive impairment, especially on memory in aging rats [56,57,58].
The thyroid hormone receptor interactor 4 (TRIP4), casein kinase 1 gamma 1 (CSNK1G1), and PCNA clamp associated factor (PCLAF; also known as PAF and PAF15) genes are located on BTA10, and they were also identified as candidate genes controlling the expression of temperament later in life [animals older than 10 years; (see Additional file 1: Table S3)]. TRIP4 gene has been identified as a candidate gene for Alzheimer’s disease susceptibility , and the CSNK1G1 gene may be a cause of syndromic developmental delay and autism spectrum disorder . Interestingly, a role for TRIP4 gene has also been identified in behavioral maturation in honeybees, which comprises the labor-division among worker bees across their lifetime . Behavioral maturation can be seen in the hierarchical emergence of activities of the work-bees over time, for example, from “2–3 weeks of adult life, workers perform tasks inside the hive, such as nursing and food storage, and as they become older, they progress to tasks outside, including foraging for pollen and nectar” .
The neuronal regeneration related protein gene (NREP; also known as C5orf13, D4S114, P311, PRO1873, PTZ17, or SEZ17) was identified as a candidate gene controlling behavioral responses at around 8 and 9 years-old (see Additional file 1: Table S3), and NREP plays an important role in the transforming growth factor beta (TGF-\(\beta\)) pathway . In pigs, the NREP gene was indicated to improve meat production based on a gene expression study comparing Czech Large White pigs and wild boars . Based on the same TGF-\(\beta\) pathway, the NREP gene has been identified as a rare variant in about 20% of suicide patients, based on a whole-exome ultra-high throughput sequencing analysis in brain samples of human suicide victims .
Genetic and genomic assessment of average temperament and learning and behavioral plasticity in cows
The Legendre orthogonal polynomial of first order was used to evaluate the average cow temperament measured at weaning and learning and behavioral plasticity, in which they are indicated by the intercept and slope of the random regression model, respectively. A similar approach was used in chickens, in which the authors showed that laying floor eggs and perching can be learned and there is a genetic component associated with it .
Learning and behavioral plasticity explains a low proportion of the total variation in cow temperament (3.33%). However, there is a genetic component associated with it, conveying that genetic and genomic selection could be applied to improve the magnitude of habituation over time, consequently, achieving the breeding goal of more docile animals. Moreover, learning and behavioral plasticity that leads to a favorable change in behavior (i.e., habituation rather than sensitization) may have a positive impact on animal welfare, both on farms that rely on frequent human–animal interaction, and high-tech farms. There is genetic variation in the learning and behavioral plasticity spectrum from habituation to sensitization among individuals, as shown in Fig. 5a. Both outcomes likely depend on well-developed cognitive abilities of learning, memory, and potentially individual human recognition, but differ in appraisal of threat. The top three sires in Fig. 5a (circle points) produced progeny that have a desirable expression of learning and behavioral plasticity because they habituate to handling. The contrasting sires (square points) produce progeny with an undesirable expression of learning and behavioral plasticity manifested by sensitization.
A wide range (from 0.02 to 0.99) of heritability estimates have been reported for cognitive performance across many species, including humans, chimpanzee, rhesus macaque, domestic pig, rat, mouse, zebra finch, honey bees, guppy, and jungle fowl [18,19,20,21,22, 61]. In this study, a low heritability was observed for the learning and behavioral plasticity component (0.02 \(\pm\) 0.002; Table 4). Strategies such as the use of genomic information [66, 67] and quantification of specific cognitive attributes (e.g., appropriate matching of threat appraisal with actual risk of harm) might assist and advance genetic improvement for learning and behavioral plasticity. Furthermore, in this study most cows had only one record, which could compromise the estimation of the slope in the RRM, consequently, rendering lower heritability for the learning and behavioral plasticity. However, we have done an additional analysis including only animals with more than two records and no major differences were observed in the estimation of the slope heritability [learning and behavioral plasticity; (see Additional file 1: Tables S10 and S11)].
The permanent environment effect explains the larger part of the variance of learning and behavioral plasticity, indicating that the experiences that the animal had throughout their life are a greater determinant of learning and behavioral plasticity between individuals. As discussed before, the environment may also induce changes in the expression (and importance) of certain genomic regions. Different rearing conditions can lead to distinct experiences, as seen in Fig. 5c–e. However, even within the same rearing conditions, animals can have individual experiences, and, consequently, different impacts on their behavioral response (shown in the animals in Fig. 5c–e).
A moderate-to-high genetic correlation was observed between the average CT and learning and behavioral plasticity. Accordingly, many of the annotated candidate genes (68%) identified overlapped between both components. The learning and behavioral plasticity had more genes uniquely identified [29 genes: e.g., (BTA1) DGKG; (BTA10) BAHD1, CCDC32, CHST14, CSNK1G1, FUT8, TRIP4; (BTA16) KISS2, PIK3C2; (BTA18) ADGR1, ADGR3, ADGR5; (BTA24) IMPA2] compared to the average CT (four genes: TOX, SMUD3, U6, 5S_rRNA).
For the average CT, two of the four uniquely identified genes were associated with mental disorders in humans and behavioral traits in livestock species. The thymocyte selection associated high mobility group box gene (TOX) has an important role in regulating neural stem cell proliferation in neural tissues , and differential expression was seen between patients with schizophrenia and control patients (accuracy of 86% ). Another gene previously described, U6, is a paralog gene associated with adrenaline levels, feeding behavior, maternal behavioral, suckling reflex, and temperament (including Angus yearling temperament) in cattle and pigs [5, 23, 24].
The cognitive-related candidate genes were previously associated with development delay, impaired learning and behavioral plasticity, and other mental and neuronal disorders in humans. The diacylglycerol kinase gamma gene (DGKG) was reported to be negatively regulated by the triiodothyronine hormone in mice inactivated for astrocyte-specific Dio2, which resulted in mood and behavioral disorders, such as anxiety-depressive behavior . Similarly, BAHD1 has also been suggested to regulate brain cells and, consequently, to be linked to anxiety-like behavior . Comparing BAHD1-knockout and control mouse groups, the authors identified many differentially expressed genes in the brain suggesting a role for BAHD1 in transcriptional regulation in neural tissue . Interestingly, 52% of the annotated cognitive-related genes are from the same family of some differentially expressed genes found in BAHD1-knockout mice, including the adhesion G protein-coupled receptor G1 (ADGRG1), coiled-coil domain-containing (CCDC family; including CCDC102A), carbohydrate sulfotransferase 11 (CHST family), casein kinase 1 gamma (CSNK1G family), and thyroid hormone receptor interactor (TRIP family) genes.
Three adhesion G protein-coupled receptor G (ADGRG) family genes had been identified as cognitive-specific candidates (see Additional file 1: Table S4), which have crucial neurodevelopment functions . Furthermore, ADGRG1 knockdown mice were associated with depression-like behavior, executive dysfunction, and poor response to neurological treatment (e.g., antidepressant ). The carbohydrate sulfotransferase 14 (CHST14), casein kinase 1 gamma (CSNK1G1), and thyroid hormone receptor interactor (TRIP4) gene families were previously associated with cognitive performance in mice and honeybees [60, 61, 74]. The CSNK1G1 and TRIP4 gene functions and associations were described above as they were implicated in CT of cows of more than 10 years of age.
Knockdown mice for the CHST14 gene were linked to spatial learning and memory impairment with the gene action located in the hippocampus . Furthermore, the CHST14 and CCDC32 genes have been identified in a genome-wide-gene-by-environment interaction as a carrier for high risk of suicide, primarily driven by post-traumatic stress in woman . The majority of the genes identified in this study sustain their function on behavior and cognitive performance across adult life. However, the CHST14 and CCDC32 genes support the hypothesis that genomic regions can change their effect and/or impact over time due to specific experiences, as shown in Fig. 8. Further validation studies on this longitudinal behavior should be carried out.
In general, the adrenal glands release two major hormones observed during stress responses: cortisol or corticosterone and aldosterone. Physiological stress-driven events activate the hypothalamic–pituitary–adrenal (HPA) axis signaling the adrenal glands. Cortisol is the primary steroid hormone to cope with stressors in some mammals . However, concomitantly with the HPA, other systems are also activated, for example, the sympathetic-adrenomedullary system leading to increased secretion of aldosterone . Excessive stressors can impair many essential mechanisms, including personality development and behavior . For instance, chronic excess of cortisol and aldosterone concentrations can lead to cognitive disorders, such as suppressed memory retrieval [77, 78]. Pathways and mediators for these steroid hormones, cortisol and aldosterone, were enriched in this study, such as cortisol metabolic, aldosterone and glucocorticoid biosynthetic, and cellular response to peptide hormone stimulus (see Additional file 1: Table S5).
Implications and next steps
The covariance structure among CT age groups was obtained from RRM fitting Legendre orthogonal polynomial of first order, in which we assumed homogeneous residual variances across age groups. This assumption was made due to computational and software limitations, but future studies should model the phenotypic variation across the years considering heterogeneous residual variance. Furthermore, the model used considered a constant permanent environment impact on the phenotype, which, for behavioral traits, may not be true. Throughout the animals’ life they could have a negative followed by a positive experience, resulting in a null permanent environmental effect because these events cancelled each other out. Therefore, it would be advantageous to evaluate correlations among ages fitting a cumulative permanent environment model for a longitudinal behavioral study . Furthermore, given that either phenotypic or genetic selection might be applied based on docility in combination with other economically important traits early in life, animals with extreme temperament might be prematurely culled and therefore, might have less phenotypic records on later age-groups (censored data). Therefore, models considering the censured nature of the data should be evaluated to account for potential pre-selection bias.
The first-order Legendre orthogonal polynomial provided the opportunity to evaluate learning and behavioral plasticity in cows over the years, which can be an indicator trait to improve the animal and handler-welfare in the long-term, as well as selection of animals that would habituate to the emerging technology on farms assuming that habituation to handling predicts habituation to other management related stimuli. Another reason for the first-order Legendre orthogonal polynomial fitted was data limitation (i.e., average repeated records per animal equal to two). Therefore, continuing to record docility scores at yearling and on older animals will be essential. Other Legendre polynomial orders could be tested in future analyses (when more records per cow are available) to identify the optimal behavioral pattern across a cow’s lifetime.
Learning and behavioral plasticity was defined by the data available; in other words, no new data recording other than already being collected by farmers was used. However, low heritability was observed providing opportunities for the identification of other optimal indicators with higher heritability within the scope of easy implementation. Furthermore, as a consequence of the high genetic correlation between overall temperament (i.e., intercept) and learning and behavioral plasticity (i.e., slope; Table 4), the current selection strategy on temperament results in indirect selection for learning and behavioral plasticity. The current genetic selection for YT would also indirectly improve the docility of cows at different ages. In addition, the inclusion of additional temperament records collected later in life into the YT genomic evaluations might increase the accuracy of the breeding values.
This is the first study evaluating the changes of genomic-region effect on cattle temperament across the animal’s life. There is variation in genomic region effect over the years for temperament (Fig. 6), and genes associated with aging-related-diseases were located in those regions. This information is key for understanding behavioral changes in a cow’s life, as well as providing opportunities to identify environmental stimulators associated with those alterations. Further longitudinal gene expression studies are key to understanding their effect on temperament, the epigenomic regions involved, and to pinpoint environmental events regulating their expression.