Introduction

Aggression in animals, including humans, serves important purposes in securing mates, territory and food and, therefore, is one of our most basic behaviors. We understand aggression as an overt behavior directed at an object or a subject with an intention to cause harm or damage (Gannon et al. 2009). This definition captures a heterogeneous and multifaceted construct, including the important distinction between two different subtypes of aggression, reactive and proactive (Smeets et al. 2017). Reactive aggression is the impulsive response to provocations, frustrations or threats from the environment. Proactive aggression is conscious and goal-directed, planned with a specific outcome in mind (e.g. attainment of resources). While both subtypes can co-occur, different neurobiological structures have been shown to underlie reactive and proactive subtypes (Blair 2013).

Animal studies have shown that reactive aggression is mediated by the limbic system, through a circuit that runs from the amygdala to the hypothalamus, and from there to the periaqueductal grey. This circuitry is regulated by frontal cortical regions, particularly by the prefrontal cortex (PFC) and the anterior cingulate cortex (Blair 2013). There is strong evidence that reactive aggression often appears earlier in life than proactive aggression, and that heightened levels of reactive aggression at early stages are associated with increased levels of anxiety in adulthood. In humans, childhood reactive aggression can be found in different neuropsychiatric conditions that can lead to conduct disorders during adulthood (Fite et al. 2010). Therefore, a comparison of gene expression profiles in PFC tissues at different stages of development may be relevant to reveal some of the molecular characteristics and mechanisms involved in aggression-related perturbations that characterize early stages of aggression and its progression into adulthood.

The Lidia cattle breed is a promising and natural animal model for studying reactive aggression, given that this breed has been selected for centuries for its agonistic, aggressive behavior by means of a series of traits registered by the breeders that classify their aggression and fighting capacity (Silva et al. 2006a, b). The breeders use as selection criteria a set of traits registered on a categorical scale, all of them showing significant heritability values that range from 0.20 to 0.36 (Silva et al. 2006a, b; Menéndez-Buxadera et al. 2017).

The Lidia breed has been raised over centuries with the exclusive purpose of taking part in socio-cultural and anachronistic events grouped under the term of “tauromaquias” (Eusebi et al. 2018). At these events, male bovines are placed in a bullring in solitary with a human subject, resulting in a series of aggressive reactions that measures the animal´s fighting capacities by means of reactive responses, such as fighting, chasing and ramming. Different types of tauromaquia take place at the bullrings, with perhaps the most popular being those known as major festivities: the “corrida” and the “novillada” (Domecq 2009). The difference between these two events is the age of the bovines. In a “corrida”, the protagonists are 4-year-old bulls, while a “novillada” uses 3-year-old bovines (Maudet 2010). This age discrimination between events is not random; 3-year-old animals are smaller, less strong and their reactive responses tend to be less explosive and intense than those of 4-year-old bovines.

Our previous studies using Lidia bovines as a model to study the neurobiological basis of aggressive behaviors revealed widespread changes in PFC gene expression compared with tamed Wagyu bovines, a meat-production breed known for their docile-temperament. Thus, we aimed to provide insights into potential molecular and cortical substrates of aggressiveness (Eusebi et al. 2020). Furthermore, genome scans for selection signatures in the Lidia cattle breed have uncovered the existence of several genes strongly associated with aggressive behaviors previously reported in humans and other animal models (Eusebi et al. 2019). Particularly, we have associated a novel polymorphism in the promoter of the monoamine oxidase A (MAOA) with different levels of aggressiveness among cattle breeds, a gene that has been clearly linked to aggression in humans (Brunner et al. 1993; Eusebi et al. 2019).

Although studies characterizing gene expression differences between juveniles and adults in the PFC have shown promising results in mice (Agoglia et al. 2017; Lander et al. 2017), to our knowledge, there is a lack of such studies using cattle as a model. Thus, the aim of the present study is to explore differences in gene expression at two age stages of Lidia cattle, “young” three-year-old and “adult” four-year-old bovines, given the marked differences in the intensity of their agonistic responses. For this purpose, gene expression profiles from the PFC tissue of eight adult and eight young Lidia bovines were compared. The results of the present study will give insights into the molecular modulation involved in the complex mechanisms leading to the maturation of the brain, with a focus on aggressive behaviors.

Material and Methods

Ethics Statement

In this study non-purposeful killing of animals was made. The samples we used belong to animals that were slaughtered within standard procedures approved by the Spanish legislation applied to abattoirs (BOE 1997) in a commercial activity. Thus, no special permits were required to conduct the research. Samples were collected from bovines after slaughter. No ethical approval was deemed necessary.

Animals

The breed´s inherent aggressiveness was considered adequate for the objectives of this study. Post-mortem PFC tissue samples were retrieved from 16 non-castrated bovine males of the Lidia breed, eight three-year-old and eight four-year-old, for the young and adult groups respectively. Animals belong to the lineages Santa Coloma-Albaserrada and Domecq, and were included in different batches according to both age and lineage (Supplementary Table 1), all affiliated to the Lidia Breeders Association (UCTL, https://torosbravos.es/). All individuals belonged to an “elite” group of aggressive bulls, selected by their breeders according to the standardized traits of aggressiveness, ferocity, face hiding and nobility on a categorical scale from 1 to 10 for each trait (Silva et al. 2006a, b). The genealogical and behavioral scores of these traits have been recorded between 1984 and 2010 and analyzed by Menéndez-Buxadera et al. (2017), using multi-trait reaction norm models, which revealed heritability values ranging between 0.230 and 0.308, with aggressiveness attaining the highest of the heritability scores. Non-related Lidia individuals were raised under an extensive farming system, pasture fed until 6–8 months prior to their sacrifice. At this stage the bulls from each batch were separated into wide-fenced enclosures and fed with a fattening supplementary diet of ad-libitum high energy and highly digestible concentrates (Lomillos and Alonso 2019).

In standard RNA-seq studies, this level of replication is 2.6-fold higher than the minimum required (3 individuals/group). Similarly to a classic resident-intruder laboratory test, where an unfamiliar animal (intruder) is placed in the territory of another animal (resident) resulting frequently in an aggressive conflict (Koolhaas et al. 2013), each bovine was placed in a bullring and incited to develop reactive aggressive responses for approximately thirty minutes prior to its sacrifice, to measure their fighting abilities (Domecq 2009). Samples of PFC were collected less than an hour post-mortem and immediately submerged in RNA-later™ (Thermo Fisher Scientific, Madrid, Spain), followed by 24 h’ storage at 5 °C and long-term conservation at −80 °C.

RNA Isolation and Sequencing

Total RNA was extracted from PFC samples using the RNeasy Lipid Tissue Mini Kit (QIAGEN, Madrid, Spain), and total RNA concentration and purity were assessed with a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Spain). For quality check, the OD 260/280 ratio was determined to be between 1.87 and 2.0. RNA integrity number (RIN), was determined using the Bioanalyzer-2011 equipment (Agilent Technologies, Santa Clara CA, USA). To guarantee their preservation, RNA samples were treated with RNAstable (Sigma-Aldrich, Spain), and shipped at ambient temperature to DNA-link Inc. sequencing laboratory (Seoul, Korea) to perform high throughput sequencing using a Novaseq 6000 sequencer (Illumina, San Diego, CA, USA). All these procedures were conducted according to the respective manufacturers´ protocols.

Individual libraries for each of the analyzed bovines (N = 16) were processed using the TruSeq Stranded mRNA Preparation kit (Illumina, San Diego, CA, USA) according to manufacturer´s instructions. All samples were sequenced in the same flowcell and lane following a pair-end protocol (2 × 100 bp). Individual reads were de-multiplexed using the CASAVA pipeline (Illumina v1.8.2), obtaining the FASTQ files used for downstream bioinformatics analysis.

Bioinformatics Analyses

Quality control and pre-processing of genomic datasets was carried out with the PRINSEQ v. 0.20.4 software (Schmieder and Edwards 2011). To improve the downstream analysis, low quality (Q < 20) and ambiguous bases (N) were first trimmed from both ends of the reads and, then, the trimmed reads were filtered by a Phred quality score (Q ≥ 20 for all bases) and read length of ≥ 68 bp. The obtained FASTQ files were mapped to the bovine reference genome (Bos taurus ARS.UCD 1.2) with the STAR Alignment v.2.7.3a software (Dobin et al. 2013), using default parameters for pair-end reads and including the Ensembl Bos taurus ARS-UCD 1.2 as reference annotation. Once reads were mapped, the expression of each mRNA and their relative abundance in fragments per kilobase of exon per million fragments mapped (FKPM) were measured with Cufflinks v.2.2.1 (Dobin et al. 2013). Transcripts with a 0 (zero) FPKM read-out in any one of the samples were eliminated due to causing division errors; only expressed transcripts that met the minimum quality values were included for further analyses, leaving a total of 14 animals (eight young and six adults).

The analysis of differential gene expression between young and adult groups was performed using Cufflinks software (Trapnell et al. 2010). The assembled transcripts from all samples were merged using the Cuffmerge command. Cuffdiff was then used applying default parameters and the option “fr-firststrand” to define pair-end reads. A Benjamini–Hochberg False Discovery Rate (FDR), which defines the significance of the Cuffdiff output, was set as threshold for statistically significant values of the Differentially Expressed Genes (DEG). The standard FDR < 0.05 was applied as cutoff. The results of the Cuffdiff RNA-seq analysis were visualized with the R software application CummeRbund v.2.28.0 (Goff et al. 2012).

To examine the relationships between differences in the PFC gene expression among groups and their biological functions, the Log2 Signal Fold Change (FC) score was used to partition the DEGs into up regulated and down regulated groups.

Cross Species Comparative Analysis (CSCA)

We conducted a comparison between our DEG and a compendium of genes associated with aggressiveness and previously identified in different genomic studies in humans, rodents, foxes, dogs and bovines, as proposed by Zhang-James et al. (Zhang-James et al. 2019). This list is based on four categories of genomic evidence: (a) genome wide association studies (GWAS) in humans of different age groups, children and adults (Fernández del Castillo and Cormand 2016); (b) genes showing selection signatures previously identified in the Lidia breed (Eusebi et al. 2018, 2019); (c) RNA-seq studies in rodents (Clinton et al. 2011; Malki et al. 2014) and silver foxes (Kukekova et al. 2011; Kukekova et al. 2018); and (d) genes identified in studies with causal evidence of the Online Mendelian Inheritance in Man (OMIM) and the Online Mendelian Inheritance in Animals (OMIA) databases, and a report of knockout (KO) studies in mice (Våge et al. 2010; Veroude et al. 2016; Zhang-James et al. 2019). The gene-list and details of the different studies are described in Supplementary Table 2. Bovine official gene names were converted to their human orthologues using biomaRt (Durinck et al. 2005) in order to homogenize our DEG with the cross-species gene list. We assigned a weight value (weighted ranking, WR) with the same conditions proposed by Zhang-James et al. (2019). For statistical analysis of the intersection between DEGs and genes identified in studies of aggression, we cross-referenced each gene list using the PANTHER (www.pantherdb.org), NCBI HomoloGene (www.ncbi.nlm.nih.gov/homologene), and Ensembl orthologue databases Bos taurus ARS-UCD 1.2 and Human reference (GRCh38.p13) genomes. If no human–bovine one-to-one orthologues were found in any database, we removed the relevant genes for statistical analysis.

Ingenuity Pathway Analysis

Up and down-regulated DEG lists were imported into the Ingenuity Pathway Analysis (IPA) (QIAGEN, www.qiagen.com/ingenuity) software to assess Gene Ontologies (GOs) and canonical pathways enrichment scores. Additionally, up and down-regulated DEGs in common with the CSCA were imported together to build biological networks and identify upstream regulators.

For the canonical pathway analyses, a right-tailed Fisher’s exact test was performed for the enrichment of the DEGs in IPA’s hand-curated canonical pathway database. Here, the P-value calculated for a pathway measures the probability of being randomly selected from all of the curated pathways. To control the error rate in the analysis, IPA also provide the corrected P-values to identify the most significant results in IPA’s canonical pathways based on the Benjamini–Hochberg method (Benjamini and Hochberg 1995). This tool allowed us to identify the signaling pathways in which the DEGs were enriched. We used the predetermined cut-off of the corrected P-value < 0.05 to define the significant pathways.

For network enrichment and detection of upstream regulators, the genes in common with the CSCA were overlaid onto a global molecular network (GMN) developed based on the Ingenuity Pathways Knowledge Base, in which functional relationships such as activation, chemical-protein interaction, expression, inhibition and regulation of binding were manually searched. Sub-networks of genes were then extracted from the GMN based on their connectivity using the algorithm developed by IPA. The sub-networks were visualized using IPA’s Path Designer tool.

Results

Sequencing and Read Assembly

For the 16 individually sequenced samples, we generated an average of 61.6 million pairs of 100-bp paired-end reads per sample after filtering by quality score. For each sample, an average of 91.6% reads were uniquely mapped to the bovine reference genome, ranging from 86.8 to 95.2% (Supplementary Table 1). The transcripts expressed in both the young and adult groups reaching FKPM values < 0 (N = 8 young and N = 6 adult) were included in the differential expression and subsequent IPA analyses (Supplementary Fig. 1).

A total of total of 27,153 genes were revealed in young and adult groups. Of those genes, only 243 were differentially expressed; 50 up-regulated and 193 down-regulated (Fig. 1A and B). For the detailed list of up and down-regulated DEGs see Supplementary Table 3.

Fig. 1
figure 1

Heatmap of up-regulated (A) and down-regulated (B) differentially expressed genes (DEG) in the adult group

Full size image

Genes in Common with the Cross-Species Comparative Analysis (CSCA)

The cross-reference analysis between our up and down-regulated DEGs in the adult group and the gene-list associated with aggressive behavior, overlapped at numerous loci. Notably, the down-regulated DEGs shared 21 genes, while the up-regulated displayed only 8 genes in common with the gene reference list (Table 1). The subset of DEGs included in the CSCA are shown in the expression bar plot of Fig. 2, which details their standard deviation of FPKM values. As well as these 29 genes, a many-to-many orthologous match was found between the bovine non-classical major histocompatibility complex class I antigen (BOLA-NC1) and the human equivalents (HLA-A, HLA-B, HLA-C, and HLA-F).

Table 1 Up and down regulated DEGs in four-year-old bulls in common with the cross-species comparative analysis (CSCA)
Full size table
Fig. 2
figure 2

Bar chart of up and down regulated DEGs in common with the cross-species comparative analysis (CSCA). Gene abundance is represented in fragments per kilobase of exon per million fragments mapped (FPKM) of the adult (4-year-old) and young (3-year-old) groups

Full size image

Statistical Analysis of Aggression-Associated Differentially Expressed Genes (DEG)

In order to test whether the 29 adult Lidia DEGs identified in the CSCA represent a statistically significant association between maturation and aggression, we calculated the cumulative hypergeometric probability of this overlap occurring. Following removal of non-orthologous genes, 1701 aggression-associated genes across species and 238 Lidia DEGs remained. Taking the estimated 22,000 genes in the bovine genome (Elsik et al. 2009), an overlap of at least 29 between these two gene samples was significantly unlikely to occur (p-value: 0.0099). We carried out the same analysis taking only aggression-associated genes that were identified in brain-expression studies across the different species (totaling 1151). Of the 238 Lidia DEGs, 21 were also identified in these studies (p-value: 0.0137). Under a more restrictive analysis, whereby only cortically expressed genes are taken as potential comparators—estimated at 85% of all protein-coding genes in both humans and pigs (Sjöstedt et al. 2020)—, the intersection of 21 was slightly more likely to have occurred by chance (p-value: 0.0617).

Ontological Analysis of the Differentially Expressed Genes (DEG)

The up and down-regulated DEGs in the 4-year-old adult group were subjected to a Fisher’s exact test within Ingenuity Package Analsis (IPA) to obtain Ingenuity ontology annotations. IPA annotations follow the GO annotation principle, but are based on a proprietary knowledge base of over a million protein–protein interactions. We used Fisher’s exact test for annotation and the FDR for multiple testing corrections, both for the up and down regulated DEG with P-values ≤ 0.05, to infer their pathway enrichment scores. Significant results are summarized in Supplementary Table 4. The most relevant results for the up-regulated DEGs were obtained under the physiological system development and the disease and disorders categories. Within the first of these categories, the top of the list gathered terms related to behavior (highest p-value range of 5.51E-08 and 15 DEG) and nervous system development and function (highest p-value range of 7.30E-07 and 24 DEG), while within the category of diseases, terms related to cancer (highest p-value range of 1.21E-07 and 44 DEGs), neurological diseases (highest p-value range of 8.07E-07 and 21 DEGs) and psychological disorders (highest p-value range of 3.35E-06 and 12 DEGs) occupied the top positions.

Additionally, we explored the overrepresentation of specific enriched functions within the category of behavior and its relationship with diseases or function annotations, finding high overall significance levels of association between up-regulated DEGs and diverse behavioral conditions (p-values < 10–5) (Table 2, Fig. 3).

Table 2 Overrepresentation of enriched specific functions within the category of behavior
Full size table
Fig. 3
figure 3

Regulator effects of the IPA package applied to the up-regulated DEG and diverse behavioral conditions (p-values < 10–5). In the lower tier, the expected behavioral consequences of the up-regulated DEG changes are shown by considering the Ingenuity Knowledge Base. In orange color are marked those functions predicted to be activated and in blue those predicted to be suppressed (Color figure online)

Full size image

The Core Analysis of the down-regulated DEGs in the 4-year-old group also revealed genes related to the physiological system development and the disease and disorders categories. The top enriched genes in the first category were associated with organismal (highest p-value range of 2.38E-17 and 63 DEGs), embryonic (highest p-value range of 2.38E-17 and 58 DEGs) and cardiovascular system (highest p-value range of 1.08E-16 and 43 DEGs) development functions. Furthermore, at the diseases and disorders category the down-regulated DEGs were shown to be involved with cardiovascular disease (highest p-value range of 3.05E-11 and 34 DEGs), cancer (highest p-value range of 1.73E-10 and 84 DEGs), and organismal injury and abnormalities (highest p-value range of 1.73E-10 and 84 DEGs) (Supplementary Table 4).

Metabolic Pathways Affected by the Up and Down-Regulated DEGs in the Adult Group

Within the up-regulated DEGs in the adult 4-year-old group, four pathways surpassed the significance threshold (p-values < 0.05, Supplementary Fig. 1). Synaptogenesis and Wnt/β- catenin signaling were among the most significant pathways (p-value = 3.08E-03). Besides, the glutamate receptor pathway, which mediates synaptic signaling by this primary excitatory neurotransmitter in the central nervous system (CNS), also resulted highly significant (p-value = 5.74E-03).

The analysis of down-regulated pathways revealed the involvement of the retrieved DEGs on 53 metabolic routes (Supplementary Table 5). Amongst the most significant of these, the intrinsic prothtrombin activation pathway (p-value = 1-08E-08) and, notably, again the Wnt/β- catenin signaling (p-values = 4.76E-05) route, are highlighted within the down-regulated DEGs in the adult group (Supplementary Fig. 1).

Gene Networks and Upstream Activator Analyses of the Genes in Common with the CSCA

We analyzed the gene transactivation networks and upstream activators of the 29 DEGs with one–to–one orthologous matches in the CSCA, finding three highly interconnected networks related to organ morphology, organismal development, behavior, cancer, neurological disease and psychological disorders (Fig. 4).

Fig. 4
figure 4

Network analysis of the 29 DEGs in the adult group in common with the CSCA. Genes highlighted in color correspond to the DEGs in our study; in red the up-regulated and in green the down-regulated DEGs (Color figure online)

Full size image

Finally, the upstream analysis tool of the IPA package was employed to identify potential upstream regulators that may explain the differential patterns of expression between young and adult groups. We explored the activation/inhibition states of the top upstream regulators by using IPA’s activation z-score tool, identifying four main upstream regulators: KRAS proto-oncogene (KRAS), Brain Derived Neurotrophic Factor (BNDF), CAMP Responsive Element Binding Protein 1 (CREB1) and SRY-Box Transcription Factor 2 (SOX2). As shown in the graphic network in Fig. 5, these genes are involved in an array of diverse biological functions related with behavioral development.

Fig. 5
figure 5

Network visualization of the four major upstream regulators of the DEGs in common with the CSCA. In orange, upstream regulators predicted to be activated; in red and green, genes whose expression increases or decreases in response to the activation of upstream regulation, respectively. The shapes of the nodes reflect the functional class of each gene product, as shown in the legend. The symbols marked in grey define functions. Solid and dashed lines between genes represent known direct and indirect interactions, respectively, with orange lines leading to activation and blue to inhibition (Color figure online)

Full size image

Discussion

The structural and molecular changes occurring in the PFC of diverse species at different age stages have been previously explored (Counotte et al. 2010; Moczulska et al. 2014;, Gonzalez-Lozano et al. 2016). The majority of these analyses conducted in rodents focus on the study of the molecular changes occurring in postnatal and puberty-adolescence stages. However, studies on the specific changes of genetic expression during adulthood are still scarce (12, Agoglia et al. 2017; Lander et al. 2017). Here, we used RNA-seq to identify genes and gene networks that show differential expression in the PFC of adult bulls when compared to young bovines of the aggressive Lidia breed.

A total of 243 significant DEGs, 50 up and 193 down-regulated in the adult group were identified. The GO analysis of up-regulated DEGs in the adult group rendered highly significant results for behavioral and nervous system development features, while the down-regulated genes were mostly implicated in systemic developmental functions (organismal, embryonic and cardiovascular system). We detected a modest coincidence between these results and what has been published in previous analysis using high-throughput molecular data in the PFC. For instance, Agoglia et al. (2017) used an unbiased proteomics approach to analyze the differential expression of proteins in the adolescent mouse PFC compared to that of adults detecting a total of 58 differentially expressed proteins. Among them, only the differential expression of the Fatty Acid Binding Protein 5 FABP5 protein/gene was also detected in our study. This gene belongs to the fatty-acid binding proteins family, whose aberrant expression has been implicated in cellular dysfunction in neurodegenerative disorders (Cheon et al. 2003). Nonetheless, there is a concordance between the proportion of up and down-regulated DEGs in the adult groups reported by both us and Agoglia et al. (2017), as well as results from the GO analyses showing down-regulation of organismal growth and development processes in adults. This suggests that both technologies, proteomics and RNA-seq are correlated, as may be expected, detecting similar sets of DE molecules.

Among the 29 genes identified in the cross-species comparative analysis (Table 2), many of them play important roles in organismal development, cellular morphology, neural plasticity and molecular transport. It is also worth highlighting two genes implicated in functional development: the NPASG gene (up-regulated DEG with prominent WR value; Table 1) participates in the structural and functional plasticity of neurons, acting as a regulator of the excitatory-inhibitory neural balance (Jaehne et al. 2015). Also, the HTR2C gene (Down-regulated DEG; Table 1) modulates neural activity, acting as down-stream effector of insulin and glucose homeostasis (Stam et al. 1994). have also been candidates in studies of human pathophysiology in mental disorders, including schizophrenia, bipolar and cognitive disorders (Iwamoto et al. 2009), and are up (NPAS4) and down-regulated (HTR2C) in the adult group of bulls.

Canonical pathway analyses also revealed that metabolic routes involved in cellular signaling and neurotransmission impact on the adult PFC during the aging process (Fig. 2). We detected an over representation of the synaptogenesis and the glutamate signaling reception pathways in the four-year-old bulls, both involved in brain’s synaptic plasticity; the synaptogenesis pathway is implicated in developmental processes of formation, maintenance and refinement of synapses (Cohen-Cory 2002), while glutamatergic signaling acts as excitatory driver of the Hypothalamus–Pituitary–Adrenal (HPA) system, a complex system that triggers reactive aggression responses (Herman et al. 2004; Evanson and Herman 2015). Similar results were obtained by Lander et al. (2017); they detected an over expression in adults of glutamate markers comparing the cortical mRNA expression of mid-adolescent and adult mice.

Most of the genes making up the glutamate receptor signaling pathway are implicated in NMDA and AMPA-mediated excitation driving long term potentiation at mature synapses. This contrasts with findings that metabotropic and kainate receptor genes—which tend to down-regulate glutamatergic excitation at synapses, including in limbic circuits that control HPA activity—are disproportionately implicated in the evolution of domesticated animals and modern human species, which display reduced reactive aggression relative to their closest extant relatives (O’Rourke and Boeckx 2020).

In our recent comparison between the Lidia and Wagyu we noted evidence for down-regulated glutamate-to-GABA conversion in the more aggressive breed, suggesting heightened excitatory signaling (Eusebi et al. 2020). In previous studies of selective sweeps distinguishing the Lidia from tamer breeds, we also noted differential signals of selection on the domestication-associated glutamate receptor gene GRIK3, further suggesting that heightened excitatory signaling may play a role in the increased reactive aggression of this breed (Eusebi et al. 2018). In future studies it may be worthwhile to compare the age-dependent PFC expression profile of tame cattle breeds with that of the Lidia. Under the hypothesis that increased excitatory signaling contributes to driving aggressive behaviors, one may expect that genes tending to contribute to a relative down-regulation of glutamatergic signaling or delayed maturation of excitatory synapses may be targeted in tamer breeds.

Interestingly, we detected differences in the expression of genes included in the Wnt//β- catenin signaling pathway in both up and down-regulated DEGs in the adult group. The complex interactions of genes, with actors triggering the expression of the Wnt//β- catenin signaling pathway while others are needed to silence routes leading to its suppression, may explain the inclusion of genes from this network in the up and down-regulated cohorts. This pathway is considered one of the most strongly conserved routes among species and is involved in synaptic transmission activities (Maguschak and Ressler 2012). Malfunction of many components of the Wnt//β- catenin signaling pathway in the adult brain can result in altered behavior and cognitive disorders (Svenningsson et al. 2002).

In addition, the prothrombrin activation pathway acts as signaling route in the CNS, inducing the decrease of cyclic AMP (CAMP) and the activation of protein kinase C (PKC), both enzymes linked to catecholamine and serotonin receptors (David et al. 2004). It has been proposed that, during encounters between individuals, a de-synchronization of the top-down control at these receptors influence a lack of control of emotional reactivity (Kretschmer et al. 2014). As suggested by Agoglia et al. (2017), the increased expression on these routes observed in the 4-year-old PFC may contribute to the fine-tuning of synaptic transmission, resulting in a superior control of the PFC executive functioning in mature brains, intensifying the addressed aggressive behaviors.

The significant overlap between the age-associated DEGs we identified and the genes associated with aggression in other species strongly suggests an important influence of maturation on the development of aggressive traits in Lidia cattle. This contrasts with the retention of a docile temperament from adolescence to adulthood in other cattle breeds and marks the Lidia as a viable model for research into the neurobiological basis of behavioral neoteny under domestication. Complementarily to these results, the IPA network tool applied to the 29 DEGs with one-to-one orthologues in the CSCA, retrieved three inter-connected networks associated with organismal morphology and development, behavior, cancer, neurological disease and psychological disorders (Fig. 3). Particularly, the Mitogen Activated Protein Kinase (MAPK) and two Extracellular Signal Regulated Kinases (ERK) appear to be involved in diverse molecular mechanisms underlying aggression. Similar results were obtained in previous studies on PFC gene expression (Zhang-James et al. 2019; Eusebi et al. 2020), identifying a differential expression of the ERK/MAPK signaling pathway in their comparisons between aggressive and non-aggressive individuals.

At the regulatory level, three of the four genes predicted to be activated display clear roles in aggressive behaviors, including BDNF (Kretschmer et al. 2014; Vigers et al. 2012), CREB (David et al. 2004) and SOX2 (Gatewood et al. 2006). The behavioral functions affected by these regulators include learning and memory cognitive processes, along with behavioral conditions and neuropsychiatric disorders, such as fear, addictive behaviors, schizophrenia and bipolar disorder, all displaying aggressiveness (Fig. 4). In this regard, there is evidence of increased memory encoding during early adulthood (Rubin et al. 1998) and also, it has been observed that adults maintain elevated plasma levels of molecules associated with abnormal aggressive behaviors (i.e. corticosterone) longer than adolescents (Romeo et al. 2006). Furthermore, some neural receptors at adult PFC have shown heightened expression levels, conferring them higher risk of addictive and impulsive behaviors (Sonntag et al. 2014). This evidence suggests that behavioral and physiological functional processes of adults differs from that of adolescents, which may reflect the ongoing maturation of the limbic structures that regulate reactive aggression responses in adulthood.

Furthermore, the presence of the KRAS gene as an upstream regulator was also observed. Although KRAS is a proto-oncogene associated principally with cancer (Zhu et al. 2008), at the transcriptional level it acts as a short-term inductor to astrocytes in response to stimulus and as a sensor that adapts cells to metabolic needs and oxidative stress in the brain (Messina et al. 2017). The reactively aggressive responses of bulls during a corrida may trigger the mechanisms implicated in oxidative stress states, which have been observed to increase with age due to a disruption in the oxidant and anti-oxidant balance (Sies 1991). Besides, in neuropsychiatric disorders such as bipolar disorder, elevated markers in blood of oxidative stress levels has been reported among adults, relating oxidative stress as mediator on neuropathological processes of adulthood (Andreazza et al. 2008).

Interestingly, whereas Agoglia et al. (2017) reported enhanced protein alterations in networks that regulate neuronal signaling, anxiety-related behavior and neurological disease in the adolescent PFC when compared to adult mice, our gene expression results found similar associations but only in adults. This different outcome may be due to differences in the age-stage comparison between studies. While in mice the biological transition of adolescence to adulthood is well studied, defining maturation stages in cattle has proved more complex, as its lifespan depends on the production purpose of each breed, among other factors (Essl 1998). Given that the average life of a Lidia sire is approximately 15 years, and the average generation sire interval pathway (sire-offspring) is 7.5 years (Fernández-Salcedo 1990), perhaps here the comparison between early and middle age adulthood would be a more correct approximation.

The present findings suggest that the cortical gene expression in 4-year-old Lidia cattle is characterized by an enhancement of molecular mechanisms involved in the increase of aggressive behaviors and neurophysiological disorders. However, despite our results suggesting cortical enhanced genetic expression associated with reactive aggression in adults, this may not be representative of other brain regions within the limbic system, also reported as important actors in this process. Future experiments to evaluate the role of these other brain regions, particularly in the context of aggressive responses, will shed further light on this age transition period.