Analysis of the brain transcriptome in lines of laying hens divergently selected for feather pecking

Background Feather pecking (FP) in laying hens reduces animal welfare and leads to economic losses for the layer industry. FP is considered a heritable condition that is influenced by dysregulation of neurotransmitter homeostasis, the gut microbiome, and the immune system. To identify genes and biological pathways responsible for FP behavior we compared the brain transcriptomes of 48 hens divergently selected for FP. In addition, we tested if high feather peckers (HFP) and low feather peckers (LFP) respond differently to light since light has been shown to trigger FP behavior. Results Of approximately 48 million reads/sample an average of 98.4% were mapped to the chicken genome (GRCg6a). We found 13,070 expressed genes in the analyzed brains of which 423 showed differential expression between HFP and LFP. Genes of uncertain function and non-coding RNAs were overrepresented among those transcripts. Functional analyses revealed the involvement of cholinergic signaling, postsynaptic activity, membrane channels, and the immune system. After the light stimulus, 28 genes were found to be differentially expressed. These included an interaction cluster of core components of the circadian clock. However, differences in the response to light between HFP and LFP were not detectable. Conclusions Genes involved in cholinergic signaling, channel activity, synaptic transmission, and immune response were found to be involved in FP behavior. We propose a model in which the gut microbiota modulates the immune system, which in turn affects cholinergic signaling. This might have an influence on monoamine signaling with possible involvement of GABA or glutamate signaling.


Background
Feather pecking (FP) is a serious problem of the layer industry causing economic losses and massive impairments of animal welfare. The propensity to perform this damaging behavior is a complex trait influenced by numerous environmental factors as well as a genetic component. Reported heritability estimates of around 0.15 [1][2][3] indicate the possibility of genetic selection against FP. Despite intensive research during the last decades, the underlying mechanisms are still not well understood.
A longstanding theory explained FP as a redirected foraging behavior [4,5], which was not confirmed in later studies. Instead, it has been linked to general locomotor activity [6,7] and feather eating [7,8]. With respect to the latter, it was reported that the inclusion of feathers in the diet affects the chicken's gut microbiome [9] as well as the actual pecking behavior [10].
Comparisons of the cecal microbiomes between layer lines divergently selected for FP or with different actual FP behavior revealed apparent differences [11][12][13]. It is difficult to determine, whether these differences are a cause for FP or just a consequence of feather consumption. The hypothesis of a link with the gut microbiome is, however, convincing as FP is influenced by the serotonergic system [14][15][16], which comprises a central and a peripheral part [17]. Central serotonin partly controls dopamine and thus also affects reward-related behaviors. The peripheral part is mainly found in enterochromaffin cells of the gastrointestinal tract. Both systems do interact via the bloodstream or the vagus nerve [17]. Gut microbiota can affect the serotonergic system in various ways such as the activation of the immune system [17,18]. The latter is notable as FP has also been directly linked to immune response [19,20]. Comprehensive genomic and transcriptomic studies might help to bring ends together by identifying potentially involved genes and pathways but are still rare. Several gene expression studies linked FP to genes involved in serotonergic signaling, either directly [21,22] or indirectly via dopaminergic [22] or GABAergic [21] signaling or via a link to obsessive-compulsive disorders [23]. Mapping studies identified loci connecting FP to dopaminergic [24], GABAergic [25] and serotonergic [2,24,26] signaling as well as to the immune system [26].
As comprehensively reviewed by [18], it thus seems likely that peripheral serotonin, the gut microbiome, and the immune system are interacting with central serotonin and dopamine to influence FP. Despite the huge body of evidence for the involvement of the mentioned mechanisms in FP, it is still not understood in detail how they interact. With the aim to gain further insight into the underlying mechanisms, we used RNAsequencing (RNA-seq) to analyze the whole brain transcriptomes of 24 full-sib pairs of hens from two layer lines divergently selected for FP behavior. As light intensity is known to trigger the actual behavior given the respective motivation [27,28], we not only compared transcriptomes between lines but also before and after light stimulation. To our knowledge, this is the first comprehensive brain transcriptome study using RNAseq with respect to feather pecking in layers. Our key results point to an involvement of cholinergic signaling and a contribution of the immune system, which complements with previous studies and might be a step towards an integrated model for FP.

Transcriptome and differential expression analyses
The average number of reads per sample was 48,502,572 (SD = 6,816,183, MIN = 36,337,624, MAX = 69,407,025) with an average mapping efficiency of 98.4% after alignment to the reference genome with TopHat (detailed summary in Additional file 1). A total of 13,070 genes were expressed in the brains of the 48 analyzed animals (average FPKM > 1, FPKM values are summarized in Additional file 2). Differential expression analysis revealed that 423 genes were differentially expressed (DE) between high feather peckers (HFP) and low feather peckers (LFP) (abs. LogFC ≥1, p adj < 0.01) (Figs. 1a and 2a). Among those, we found 305 genes of uncertain function (LOC symbols), of which 156 are classified as non-coding RNAs (ncRNAs). Light stimulation led to the discovery of 28 DE genes (abs. LogFC ≥0.5, p adj < 0.01) (Figs. 1b and 2b) with no differences between HFP and LFP (summary of differential expression analysis in Additional file 3). The top 40 DE genes with the lowest p adj between HFP and LFP are plotted in a Heatmap (Fig. 3 Detailed information on putative eRNAs are available in in Additional file 4.

Functional and pathway analyses
The results of the Gene Ontology (GO) and KEGG pathway analyses with DE genes between HFP and LFP hens with clusterProfiler [29] are shown in Table 2 (complete results in Additional file 5). Most of the GO terms assigned to the comparison of HFP and LFP lines are linked to synaptic functions and membrane channels when DE genes with an absolute LogFC < 0.5 were neglected. In particular, the DE genes CHRNA2, CHRN A9, CHRNB3, and CHRNB4 were assigned to these GO terms.
Concerning the comparison of DE genes between animals before and after light stimulation (Table 3), the majority of IDs belonging to the category GO biological processes were linked to circadian rhythm. The DE genes ARNTL, NFIL3, NPAS2, and PER2 are overrepresented in these GO terms. The category GO molecular function exclusively contains IDs related to gene expression.
Geneset analysis with STRING [30] was performed on DE genes between HFP and LFP lines, where genes with an absolute LogFC < 1 were neglected (Table 4, complete results see Additional file 6). The majority of terms found in this analysis are linked to the immune system and immune response.
Analysis of DE genes before and after light stimulation with STRING (Table 5) led to the discovery of terms related to biological/circadian rhythm and gene expression as well.
Protein interaction network analysis with STRING revealed one large and one smaller interaction cluster in the comparison of HFP and HFP DE genes (Fig. 4a). The large cluster, consisting of the DE gene FIP1L1 as well as SSU72, CPSF2, CPSF3, CPSF4, WDR33, CDC5L, PCF11, CLP1, CSTF1, CSTF2, CSTF3, and PAPOLA, exclusively contains genes, which are responsible for mRNA processing as well as mRNA binding. The slightly smaller cluster with fewer connections exclusively contains genes, which are associated with the immune system -T-cell proliferation and regulation in particular. The immune cluster consists of the DE genes CD28 and LIF as well as CD274, CTLA4, CD80, CD86, IL2B, IL12A, IL6ST, and LIFR. Analysis of DE genes before and after light stimulation led to the discovery of one interaction cluster, which contains the DE genes ARNT L, PER2, PER3, NFIL3, NPAS2, and LONRF3, all of which, except LONFRF3, are core components of the circadian clock (Fig. 4b).

Discussion
In this study, we compared 24 full-sib pairs of laying hens selected for HFP and LFP in a transcriptomics approach. Besides the line comparison, we assessed the transcriptional changes before and after light stimulation. To our knowledge, this is the first brain transcriptomic study in FP that used RNA-seq. Comparison with the few other available transcriptome studies revealed only little overlap, which might be due to the fact that these studies are microarray-based and applied different designs. Hughes and Buitenhuis [31] used a more complex phenotyping scheme distinguishing between severe and gentle FP, Brunberg et al. [23] particularly analyzed hypothalamic RNA and Wysocki et al. [22] considered much older birds. Some genes, however, that these authors reported were also found DE in the current study. Notably, Hughes and Buitenhuis reported the expression of IRF2 (interferon regulatory factor 2) to be associated with severe FP This gene plays a major role in immune response [32] and might match our results with respect to immune related genes. In our study, however, we did not identify this gene, but found IRF5 (p adj = 0.0002), IRF6 (p adj = 1.17E-06), and IRF8 (p adj = 0.0002) to be DE between the lines. Wysocki et al. reported DE of genes related to monoamine signaling, namely MAOA (monoamine oxidase A) and HTRR1B (5-hydroxytryptamine (serotonin) receptor 1B). We were able to confirm MAOA (p adj = 7.2E-05), but not HTRR1B. Instead, we found DE of HTR1A (p adj = 0.04), HTR1E (p adj = 0.001) and HTR6 (p adj = 2.34E-05).
Most of the GO terms found enriched in the current study in genes DE between HFP and LFP lines are linked to synaptic function when a LogFC threshold of 0.5 is used. This, in particular, is due to the DE genes CHRN A2, CHRNA9, CHRNB3, and CHRNB4 coding for subunits of nicotinic acetylcholine receptors (nAchR). Most other studies found serotonergic or dopaminergic signaling to be involved in FP [18] and we recently re-  based on genomic studies in the very same population as used in the current study [21,25]. Cholinergic signaling has not been implicated in FP so far, but is linked to both monoamine as well as GABAergic signaling and is implicated in certain psychiatric disorders [33]. In particular, α7 nAChRs have a crucial role in the dysfunction of cortical parvalbumin-positive GABAergic neurons, as seen in schizophrenia [34]. The activation of the dopamine system by stimulating α7 nAChRs is crucial in drug dependence and withdrawal effects including motor function and rewarding properties of drug intake [35,36]. In turn, it was shown that α7 nAChR blockade in mice can reduce anxiety-like behaviors as a consequence of nicotine withdrawal [37]. Nicotinic α4β2 receptors, on the other hand, have been implicated in obsessivecompulsive disorder and it has been shown that positive allosteric modulation of these receptors attenuated the obsessive behavior [38]. The actual behavioral effects of nAChRs are mediated by dopamine and glutamate and ultimately by GABAergic neurons [39,40]. Very recently, it was shown in rats that β2 nAChRs on dopamine neurons in the ventral tegmental area mediate nicotine's conditioned aversive effects, while those on GABA neurons mediate the conditioned rewarding effects [41]. In our study, we did not find the same subunits as described above to be differentially expressed, but all other studies so far have been conducted in humans, mice or rats and little is known in chicken.   be on top of this signaling hierarchy. In a second gene set analysis, we used a more stringent set of genes by setting a LogFC threshold of 1.0. Interestingly this led to a shift in enriched terms compared to the aforementioned analysis to immunoglobulin-associated genes. We detected a protein interaction cluster consisting of 10 genes (the DE genes CD28 and LIF as well as CD274, CTLA4, CD80, CD86, IL2B, IL12A, IL6ST, and LIFR), all of which are linked to T-cell proliferation and regulation (Fig. 4a). In a recent study on post-mortem human brains, an increased number of T-cells was found in the brains of schizophrenia patients [42]. A probable  mechanism, by which an increase in the number of T cells could influence the propensity to FP has been discovered in a murine cell culture system. Mashimo et al. demonstrated that T cells are able to synthesize acetylcholine (ACh) but also respond to it via nAChRs by increased cell proliferation and increased Ca 2+ levels [43]. We propose a mechanism in which an overrepresentation of T cells leads to excess ACh production in the brain that has two major effects: (i) nAChR activation in neurons and (ii) further increase of T cell proliferation in a positive feedback loop. Since the immune-modulating genes show a higher LogFC compared to genes encoding nAChR subunits, we hypothesize that the dysregulation of T cell proliferation and activation is the initial cause of this neuropsychiatric phenotype. A recent extensive transcriptome study in humans suffering from autism spectrum disorder, schizophrenia, and bipolar disorder identified pathways related to the immune system and to transmembrane transporters/receptors to be differentially regulated in comparison to healthy individuals [44], which indicates that similar mechanisms are responsible for neuropsychiatric disorders between species. Interestingly these researchers found a large amount of DE ncRNAs between patients and controls, which they termed "psychiatric ncRNAs". With 36.9%, ncRNAs make up a large part of significant DE transcripts between HFP and LFP hens. Gandal et al. already proposed these putative "psychiatric ncRNAs" to be responsible for transcriptome dysregulation by regulating local splicing events. This might be connected to the observation that HFP animals show a globally reduced variance in gene expression [31]. However, comparison of cDNA sequences of the DE ncRNAs identified in this study with the "psychiatric ncRNAs" from Gandal et al. with discontiguous megablast [45] showed sequence homology between only one pair of ncRNAs: The human ncRNA ENSG00000234773, which is a novel zinc finger protein pseudogene, is partially homologous to LOC112533587. Due to the taxonomic distance, the huge difference in genome size, and the fact that ncRNAs are less evolutionary conserved than protein coding genes [46], we did not expect to find homologous ncRNAs between Homo sapiens and Gallus gallus. By finding the closest DE expressed protein coding genes to DE ncRNAs between HFP and LFP chickens we were able to identify 21 putative eRNAs (Table  1). Among possible targets of these eRNA candidates are ARHGEF38 and SH3YL1, which have been linked to bipolar disorder and suicide [47,48], GLRA1 and MCTP2 were found to be associated with schizophrenia [49,50], and CHIR-B3, MHCIA7, MHCIY as well as MICA are genes involved in the immune system [51]. Thus, further experiments, like overexpression studies, will help to  [28]. Light stimulation led to DE of five genes, which are core components of circadian clock (Fig. 4b) but there were no differences detectable in the light response between HFP and LFP. Whereas in the current study whole brains where analyzed, it could be especially useful to restrict analyses to particular brain regions e.g. arcopallium or caudocentral nidopallium where differences in the serotonin or dopamine metabolism between LFP and HFP were found [52,53]. Furthermore, analyses should be conducted in a larger number of birds and accounting for the actual phenotype, i.e. comparing performers against nonperformer. Here, we compared divergently selected lines, which might also differ in other traits than FP [54]. Although the lines were initially derived from the same founder population and selection was solely based on the estimated breeding value for feather pecking, genetic drift may have caused allele frequency differences between the lines [25]. Here, actual behavior was highly correlated with line with on average 1.6 bouts per bird in LFP and 12.7 bouts per bird in HFP animals during the standardized observation time. A phenotype-based approach would have resulted in an almost identical group assignment. Futures studies should aim to validate the results in a purely phenotype-based approach using a single and unselected field population. This is, however, not trivial, because the occurrence of FP cannot be reliably predicted in commercial flocks and phenotyping is challenging under field conditions, which is the reason for many studies being conducted in selection lines.

Conclusions
The analysis of whole-brain transcriptomes in 24 full-sib pairs of hens from two White Leghorn layer strains divergently selected for FP behavior revealed differentially expressed genes and pathways related to cholinergic signaling and immune response. Thus, we hypothesize an involvement of these genes and pathways in the development of FP. Previous studies unambiguously identified monoamine signaling as the key component in FP behavior, but also pointed to an involvement of GABAergic signaling, immune response and the gut microbiota. Our results provide a first building block for an integrated model of FP. In this model, the gut microbiota might be involved in the interaction between the peripheral and the central serotonergic system as well as in the modulation of the immune response. Via the latter route, there could be a connection with cholinergic signaling, which in turn projects into monoamine signaling, directly or via GABA or glutamate signaling. In coming studies, the current results have to be validated using alternative approaches and appropriate functional studies need to be designed and conducted to further challenge and develop this model.

Experimental animals and sample collection
Experiments were conducted in White Leghorn layer strains divergently selected for FP behavior [2,55]. Strains were developed from a single founder population and selection over 15 generations was solely based on estimated breeding values for feather pecking [17]. These lines were created and are maintained at the Hohenheim University and neither commercially obtained nor from a private source. Details on hatching, rearing, and husbandry have previously been described in detail [25,56]. A total of 48 hens comprising 12 full-sib pairs from each strain (high feather peckers, HFP vs. low feather peckers, LFP) were phenotyped according to established protocols at 27 weeks of age [25]. Briefly, animals were observed on four consecutive days in sessions of 20 min by at least six different experienced observers. FP was defined as non-aggressive severe pecks or pulls directed to the plumage of conspecifics [57]. A series of pecks delivered in a short sequence without changing the behavior was recorded as a single occurrence called a bout per bird. Subsequently, birds were kept under low light conditions to prevent FP (Fig. 5). One bird from each full-sib pair was then sacrificed and brains were immediately collected for RNA isolation. Chickens were CO 2stunned and sacrificed by ventral neck cutting. The remaining birds were kept under increased light intensity (≥100 lx) for several hours until they clearly showed FP and were then sacrificed as well and brains were collected for RNA isolation. All brain samples were kept in RNA later and frozen at − 80°C until further processing.

RNA isolation
Whole chicken brains were pulverized with a Retsch mixer mill MM400 in liquid nitrogen at a frequency of 30 Hz for 1 min and immediately frozen at − 80°C. For RNA isolation 900 μl Qiazol (QIAGEN), 5 μl DX reagent and 20-30 ceramic beads were added to10-30 mg of pulverized brain. The mixture was placed in a Bead Ruptor (Omin Inc.) and homogenized with the following program: 4.5 m/sec, 3 times 15 s, 10s dwell. After 5 min incubation at room temperature 100 μl gDNA Eliminator solution (QIAGEN) was added and the sample was vortexed for 15 s. After addition of 180 μl chloroform the sample was vortexed for 15 s and incubated at room temperature for 3 min. Sample was centrifuged at 12, 000 g and 4°C for 20 min and aqueous phase was transferred to a fresh reaction tube. An equal amount of 70% ethanol was added to the sample and vortexed for 15 s. RNA was purified using the RNeasy Plus Universal Mini Kit (QIAGEN) according to the manufacturer's protocol and stored at − 80°C.

NGS library preparation and sequencing
Input total RNA was visualized on a TapeStation 4200 (Agilent) and all but three samples showed excellent RNA quality with RIN-scores > 8 (three samples with RIN scores of 6.8, 6.9 and 7.7 respectively). Samples were quantified on a Qubit 2.0 Fluorometer (ThermoFischer) and 500 ng of total RNA was used as input for the TruSeq stranded mRNA library kit (Illumina) following the manufacturers manual. Resulting libraries showed a fragment size distribution of around 300 bp and were 2 × 75 bp paired-end. Sequencing was conducted on the HiSeq 4000 (Illumina) with 10 samples / lane.

Transcriptome analyses
Raw sequencing reads were processed with Trimmomatic version 0.36 [58] for adapter removal, trimming of low quality base calls, and removal of low-quality reads. Trimmomatic was used with the settings: PE -phred33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MI NLEN:36. Read pairs were discarded if one read did not survive quality control. Trimmed reads were aligned to chicken genome version GCF_000002315.5 (see Availability of data and materials) using TopHat version 2.1.0 [59] with the settings --no-novel-juncs --min-isoformfraction 0.0 --min-anchor-length 3 -r 192 and GCF_ 000002315.5.gff (see Availability of data and materials) as known transcript file. Genomic features were extracted from the general feature format file and grouped with the R package GenomicFeatures (Version 1.40.0) [60]. summarizeOverlaps from the R package Genomi-cAlignments (Version 1.24.0) was used to count exon spanning reads [60]. Differential expression between experimental groups was analyzed with DESeq2 (Version 1.28.1) [61]. Volcano plots of differential expression analyses were created with the R package EnhancedVolcano (Kevin Blighe, Sharmila Rana and Myles Lewis (2019). EnhancedVolcano (Version 1.6.0). All R packages were obtained with Bioconductor version 3.11.

Identification of eRNA candidates
One BED file was created containing DE ncRNAs and one BED file containing all protein coding genes annotated in chicken genome version GCF_000002315.5 (see Availability of data and materials). The closest coding gene to each respective ncRNA was identified with BEDtools (Version 2.27.1) [62].

Functional analyses
Gene set analyses were conducted using the R package clusterProfiler (Version 3.16.0) [29] along with the chicken genome annotation org.Gg.eg.db (Version 3.11.4) [63] and using default settings in terms of gene set size and p-value correction. Genes differentially expressed with an absolute LogFC ≥0.5 were analyzed for enrichment in GO terms and KEGG pathways against the background of all expressed genes in our data set. For the gene set analysis with STRING (Version 11.0) [8], DE genes between HFP and LFP laying lines with an absolute LogFC ≥1 and an adjusted p-value below 0.01 were used. Protein interaction networks for these genes were also computed with STIRNG using the following settings: evidence-based coloring of network edges, inclusion of all interaction sources, medium required interaction score 0.4, max number of interactions to show 1st shell: no more than 10 interactions, 2nd shell: no more than 10 interactions. For DE genes before and after light stimulation LogFC threshold was 0.5. Settings for protein interaction networks were as follows: evidence based coloring of network edges, inclusion of all interaction sources, medium required interaction score 0.4, max number of interactions to show 1st shell: none / query proteins only, 2nd shell: none. The DFG Competence Centre for Genome Analysis Kiel (CCGA) is greatly acknowledged for the cooperation in RNA-sequencing. We further acknowledge support by the Open Access Publication Funds of the Göttingen University. We thank Judith Beier and Hanna Iffland for participation in animal slaughtering and sample preparation. The data presented here has been previously presented at the ISAG 2019 conference.
Authors' contributions CFG performed all bioinformatics analyses and wrote the manuscript. AM participated in animal slaughtering and brain extraction, and performed the RNA extraction. SP participated in animal slaughtering and revised the manuscript. SF performed the sequencing. WB and JB developed the project outline. JT developed the project outline, developed data analysis strategies, and wrote the manuscript. All authors have read and approved the manuscript.

Funding
The study was funded by the German Research Foundation (DFG) under file numbers TE622/4-2 and BE3703/8-2. The funders had no role in study design, data collection and analysis, and interpretation of data and in writing the manuscript. Publication fee was covered by the Open Access Publication Funds of the Göttingen University. Open Access funding provided by Projekt DEAL.

Ethics approval and consent to participate
The research protocol was approved by the German Ethical Commission of Animal Welfare of the Provincial Government of Baden-Wuerttemberg, Germany (code: HOH 35/15 PG, date of approval: April 25, 2017).

Consent for publication
Not applicable.