RNA sequencing demonstrates large-scale temporal dysregulation of gene expression in stimulated macrophages derived from MHC-defined chicken haplotypes.

Discovering genetic biomarkers associated with disease resistance and enhanced immunity is critical to developing advanced strategies for controlling viral and bacterial infections in different species. Macrophages, important cells of innate immunity, are directly involved in cellular interactions with pathogens, the release of cytokines activating other immune cells and antigen presentation to cells of the adaptive immune response. IFNγ is a potent activator of macrophages and increased production has been associated with disease resistance in several species. This study characterizes the molecular basis for dramatically different nitric oxide production and immune function between the B2 and the B19 haplotype chicken macrophages.A large-scale RNA sequencing approach was employed to sequence the RNA of purified macrophages from each haplotype group (B2 vs. B19) during differentiation and after stimulation. Our results demonstrate that a large number of genes exhibit divergent expression between B2 and B19 haplotype cells both prior and after stimulation. These differences in gene expression appear to be regulated by complex epigenetic mechanisms that need further investigation.

Introduction

Discovering genetic biomarkers associated with disease resistance and enhanced immunity is critical to developing advanced strategies for controlling viral and bacterial infections in various species.

Disease resistance and susceptibility depends on a variety of factors including genetics. In numerous species, disease resistance has been associated with major histocompatibility complex (MHC) haplotype, as well as polymorphisms in several immune genes such as TGFβ and TNFα[1,2]. Cytokine production, specifically secretion of pro-inflammatory molecules, has also been associated with increased resistance against disease [3,4].

Studies have demonstrated association of MHC-B haplotype in chickens and resistance to a variety of viral pathogens, including AIV, Marek’s disease virus (MDV), avian leukosis virus, Newcastle disease virus and Rous sarcoma virus [5-10] as well as other pathogens [11,12]. B2 haplotype chickens are more resistant to avian coronavirus infection than B19 haplotypes and these differences in disease resistance were observed early after infection in our previous studies [10]. This suggests that innate immunity plays a major role with the macrophage being a key player in this enhanced immune response as evidenced by the B2 haplotype birds’ greater capability to produce nitric oxide (NO) in response to IFNγ and Poly I:C [13]. In addition, B2 macrophages activated T cells more efficiently than macrophages derived from B19 haplotypes [14].

Macrophages are directly involved in cellular interactions with pathogens and demonstrate distinct immune responses from more disease resistant animals in response to infection [15-20]. In addition, macrophages release cytokines activating other immune cells and antigen presentation to cells of the adaptive immune response [21-23]. It has become increasingly clear that dysregulation of macrophage function is involved in inflammatory disease processes such as rheumatoid arthritis, inflammatory bowel disease and cancer [24-26]. Involved in these interactions are crucial molecules such as Toll-like receptors (TLRs) that recognize invading microorganisms, resulting in communication with the adaptive immune system such as increased expression of MHC surface molecules, T cell receptors and secreted cytokines [21,23]. Genetic differences in any of those molecules can potentially account for differences in immune competence and thus provide potential immunogenetic markers for disease resistance to various pathogens.

IFNγ is a potent activator of macrophages and increased production has been associated with disease resistance in multiple species [27-31]. These findings indicate that chickens with enhanced IFNγ production are more resistant to certain infections. IFNγ enhances macrophage activation, expression of MHC and nitric oxide release which aides in killing of pathogens and also increases activity of cytotoxic T cells and secretion of Th1 cytokines [31,13], underscoring how crucial this process is for innate immune competence.

Macrophage TLRs appear to be primed by IFNγ, reprogramming cellular responses to other cytokines, such as type I interferons and IL-10 and activating the Jak-STAT pathway (Janus kinase and signal transduction and activator of transcription) [24, 32, 33]. IFNγ, which increases TLR receptor availability for interaction with its ligands, has been shown to induce TLR2, 4, 6 and 9 [34-37].

The response of macrophages to an immune stimulus is not just dependent on cell surface receptor and cytokine expression. Other factors include the differentiation of monocytes into functional macrophages, a tightly regulated process that influences immune competence [38]. Recent studies demonstrated a critical role for molecules such as A2B adenosine receptor for differentiation and proliferation of monocytes and macrophage function in immunity and inflammation [39, 40]. A2B expression is induced by IFNγ and leads to increase of anti-inflammatory signaling counteracting the inflammatory response activated within the IFNγ pathway.

Taken together, these studies emphasize the genetic basis of the activation of macrophages by IFNγ playing an important role in the innate immune response signaling and providing resistance to disease. In addition to inflammatory signaling, a number of transcription factor pathways and epigenetic mechanisms all contribute to immune function. Dysregulation of any of these events will lead to an impaired innate immune response and consequently, increased susceptibility to disease.

Using an ex vivo model, we investigated the gene expression in macrophages from haplotypes B2 and B19 during differentiation and after stimulation with IFNγ. Our experimental design leveraged an initial 6 day window for monocytes to differentiate into macrophages, which was followed by IFNγ stimulation between 1 and 24 h to further characterize subsequent RNA gene expression and the molecular basis for dramatically different nitric oxide production and immune function between the B2 and the B19 haplotype chicken macrophages

Material and methods

Experimental animals

Animal protocols were performed under the approval of the Institutional Animal Care and Use Committee at Western University of Health Sciences, Pomona, California (WesternU). Fertilized eggs, descended from Modified Wisconsin Line 3, were obtained from Dr. W. Elwood Briles, Northern Illinois University, and incubated and hatched under standard conditions at (38°C/50-65% humidity) [10,13] at WesternU. In addition to daily health monitoring, fresh food and water were provided ad libitum. Experimental animals were euthanized by insufflation of isoflurane gas (Butler, Dublin, OH).

Peripheral blood collection

Whole blood samples were collected via jugular venipuncture in EDTA tubes from age matched chicks at 14–18 weeks old. At no time did the amount of blood harvested from each animal exceed 1% of body weight.

Peripheral blood mononuclear cell (PBMC) isolation

Peripheral blood mononuclear cells (PBMCs) from individual birds were isolated using the differential centrifugation as previously described [41, 42, 13] with slight modifications. Briefly, blood was mixed with an equal volume of phosphate buffered saline (PBS) and slowly layered 2:1 on a Ficoll-Hypaque gradient (density 1.083) (Sigma-Aldrich, St. Louis, MO). Samples were centrifuged for 35 min (400 x g; 23°C; brake off) for retrieval of mononuclear cells. Isolated cells were washed 3x in 10 ml PBS at low speed to remove thrombocytes (180 x g; 10 min, 23°C), counted and viability confirmed based on the exclusion of 0.1% trypan blue dye (≥ 90%). PBMCs were re-suspended in PBS to a final concentration of 5 x 107 cells/ml.

Macrophage cell culture

One milliliter of PBMC suspension (5 x 107 cells/ml) was incubated (37°C/5% CO2) for 3 h in each well of a 12-well plate containing RPMI w/o Phenol Red supplemented with, 10% heat inactivated fetal bovine serum (FBS); non-essential amino acids, (0.1mM/ml) (Invitrogen, Carlsbad, CA), L-glutamine (2 mM) (Sigma-Aldrich, St. Louis, MO), 2-mercaptoethanol (55 μM/ml) (Sigma-Aldrich, St. Louis, MO), penicillin (50 U/ml) (Invitrogen, Carlsbad, CA), and streptomycin (50 μg/ml) (Invitrogen, Carlsbad, CA). Following removal of non-adherent cells with warm PBS, medium was replenished and cells were incubated for differentiation with the exception of the -6 day sample which was lysed with 400 μl Trizol (Thermo Scientific, Waltham, MA) and stored at -80°C. Prior to the replacement of medium, adherent cell cultures were washed in warm PBS. Monocytes were cultured for 6 days to allow maturation and differentiation of cells; with medium changes occurring every 3–4 days thus ensuring that optimal nutrient requirements were met. Additionally, -3 day (t-3) samples were lysed with Trizol and stored at -80°C. The morphology of adherent cells was observed daily under bright field microscopy (20x objective).

Purity of monocyte cultures using this culture method was confirmed by IFA and FACS using monoclonal antibody KUL01 as previously described as part of a different aspect of this study [13].

IFNγ stimulation

A 50 ρg/ml ch-IFNγ solution (Invitrogen, Carlsbad, CA) was prepared in RPMI w/o Phenol Red culture medium (Invitrogen, Carlsbad, CA). After washing the cells twice with warm PBS, macrophage cultures were stimulated with 1 ml of RPMI-ch-IFNγ mixture [13].

Nitric oxide assays

Nitric oxide production was measured [10, 43, 44] to confirm macrophage stimulation in assays by interferon (data not shown). Stimulation was evaluated as yes/no based on previously published results from B2 and B19 IFNγ stimulated macrophages (10)

Sample collection and RNA sequencing

A total of 145 gigabytes of RNA sequence data was obtained from B19 and B2 haplotypes of chickens. Two birds from each haplotype were selected for inclusion in the sequencing. Each bird provided blood for extraction and isolation of peripheral blood mononuclear cells. Purified monocytes were cultured for differentiation and cell samples were collected from nine time points for each bird. Samples were collected for sequencing on the day they were cultured (Day t-6), as well as on Day -3 (t-3), Day 0 (called 0 hours), and then six additional times over a 24 hour period corresponding to 1 hour, 2 hours, 4 hours, 8 hours, 16 hours and 24 hours after interferon stimulation. Cells were lysed in wells with RLT buffer containing beta-mercaptoethanol (Qiagen, Valencia, CA) and stored at -80°C. RNA was processed with the Qiashredder and RNAeasy kit from Qiagen (Valencia, CA) according to manufacturer’s instructions and sent on dry ice to Dr. Calvin Keeler at the University of Delaware for generation of libraries and sequencing with an Illumina HiSeq 2000.

An RNA sequence library was constructed from purified RNA. The library was fragmented in order to generate appropriately sized RNA fragments suitable for templates in random primed first-strand cDNA synthesis. Second strand synthesis was completed in accordance with specifications for sequencing with Illumina’s HiSeq2000 platform.

The samples corresponding to each time point from each bird were sequenced and the data was stored in a unique file for each sequenced sample and time point. Forty FASQ files were generated from the data totaling 145 gigabytes. The average file size was 3.65 gigabytes and the standard deviation was 2.25 gigabytes. The sequencing data provided 933,107,885 reads across the biological samples and time points (Table 1). Across all time points for the two B2 samples, one produced 298,903,517 reads and the other produced 165,589,594 reads. Similarly, across all time points, the B19 samples produced 285,392,384 reads and 183,222,390. For each time point (across all four birds) sequencing reads ranged from a low of approximately 78 million reads to a high of just over 171 million reads, with most time points producing over 88.4 million reads each and a few producing over 100 million reads each.

Table 1

Sequencing reads across biological samples and time points.

		t-6 days	t-3 days	t0 h	t1 h	t2 h	t4 h	t8 h	t16 h	t24 h	TOTAL t-6d_to_t24h
BIRD-A(B2)	reads	25470757	13671735	8242127	23744417	47908759	78741373	71211600	23314834	6597915	298903517
	aligned reads	11471896	8872332	5508010	16979440	35666816	56923920	38573825	10060460	3421906	187478605
	multiple alignments	148334	43889	21791	70270	131735	229889	141032	49269	16332	852541
BIRD-B(B2)	reads	1459795	28536965	24829151	29604495	1806807	5394361	25715126	17422945	30819949	165589594
	aligned reads	87027	21502194	17672812	22534847	1268878	3911819	18816906	13460540	23620025	122875048
	multiple alignments	425	92783	63653	78722	6541	18257	86313	39605	92277	478576
BIRD-C(B19)	reads	44587849	32008868	21166414	25722656	23055824	25946021	59435090	29226936	24242726	285392384
	aligned reads	17419206	12412811	14534824	16720019	15511707	17668026	41819414	16704419	16847457	169637883
	multiple alignments	177899	114518	78430	99577	73310	65676	177083	75290	63288	925071
BIRD-E(B19)	reads	28452005	6000451	23444741	20988785	22333873	21392153	14715694	18485412	27409276	183222390
	aligned reads	10587715	775587	14000961	14095140	15379919	15685618	10520441	12590915	20507446	114143742
	multiple alignments	212095	1803	80898	73748	82985	28868	39248	62596	61264	643505
ALL 4 BIRDS	reads	99970406	80218019	77682433	100060353	95105263	1.31E+08	171077510	88450127	89069866	933107885
	aligned reads	39565844	43562924	51716607	70329446	67827320	94189383	109730586	52816334	64396834	594135278
	multiple alignments	538753	252993	244772	322317	294571	342690	443676	226760	233161	2899693

Mapping reads to reference genome and identification of splice junctions / exons

The chicken reference genome WASHUC2, corresponding to Ensembl release 70, was downloaded from Ensembl.org (http://www.ensembl.org/info/data/ftp/index.html). Annotation files included the small RNA annotation files obtained from miBase release 19 (http://www.mirbase.org/). Sequenced reads were filtered to remove low quality sequences from the data. Filtered sequences were aligned to the reference genome using Bowtie and Tophat, available along with the software package Cufflinks, from John Hopkins University Center for Computational Biology (https://ccb.jhu.edu/software.shtml). The aligned reads generated by Bowtie produced gapped alignments on the reference genome which Tophat used to identify splice junctions flanking exons. The resulting aligned reads were analyzed by Cufflinks to construct transcripts corresponding to mRNA sequences. Next, Cufflinks was employed to estimate transcript specific expression levels across the transcripts and genes within the reference genome based on the number of sequence reads for each mRNA. The sequence read data was normalized using the fragments per kilobase of transcript per million mapped reads (FPKM) method to more accurately determine expression levels. The resulting transcriptome data was loaded into the MySQL relational database to more effectively manage, explore, mine and annotate the data.

Hierarchical clustering of genes and production of heat map visualization

Gene expression data was hierarchically clustered using 1-Pearson correlation on the rows and keeping the column order conserved. The resulting clustered data set was visualized as a heat map with black representing lack of gene expression, and darker shades of blue indicating lower expression values. Dark purple represents higher expression values than any shade of blue while bright pink represents the highest expression values. For visualization purposes, the heat maps were generated with maximum heat map color assigned to expression set lower than the absolute maximum expression value contained in the entire data set, subsequently all values of expression greater than or equal to the assigned expression threshold (for example, 1000) shared the same color on the heat map (regardless of whether the actual expression level was 1000, 2000, 20,000 or 90,000). This setting provided the optimal visualization of both high and low expressed genes in the heat maps.

Gene enrichment calculations were performed using the DAVID bioinformatics database tool version 6.8 (https://david.ncifcrf.gov/). The analysis was performed using comparisons of successive time points within the B2 haplotype data to identify sets of genes that were enriched (p-value < 0.05). The B2 haplotype represents the robust macrophage phenotype as characterized by nitric oxide production compared to the B19 haplotype. Subsequently, the gene enrichment was performed on the B2 data. Gene enrichment was determined using three distinct databases: gene ontology biological process, KEGG pathways, and reactome pathways corresponding to S2, S3 and S4 Tables respectively. Because a large number of enrichment annotation terms were produced, a subset of representative highlights from each of these three enrichment analyses was chosen for inclusion in the results. Highlights were selected to provide examples of the biological process annotations, KEGG pathways annotations and reactome annotations.

PCR validation of target genes

Realtime PCR was performed on a selected number of target genes to validate RNA sequencing results. RNA was taken from macrophages stimulated with IFNγ as described above for 2 and 4 hours, unstimulated samples (0h) served as control. For Realtime RT-PCR, cDNA synthesis was performed using SuperScript III First Strand Synthesis kit (Invitrogen, Waltham, MA), according to manufacturer’s instructions. PCR conditions were as follows: 95°C for 10 min-hot start, 40 cycles of 95°C for 15 sec, 60 or 63°C depending on gene (see primers) for 30 sec according to manufacturer instructions for the Biotool 2x Sybr Green qPCR Mix (Biotool, Houston, Tx). Primer sequences were designed using Primer 3 (ATP6VOC, LITAF, IL18R, TLN-1.) Primer sequences previously published were used for TLR2, TLR4, TLR5, TLR6 and TLR7 [45]. ATP6VOC (annealing 60°C) forward TGTTGTCATGGCGGGTATTA, reverse ACAAATAACCTGGGCTGCTG; LITAF(annealing 60°C) forward ATCCTCACCCCTACCCTGTC, reverse GACGTGTCACGATCATCTGG; IL18R (annealing 63°C) forward CTCTTCGTGCCTCCATTGAT, reverse ACCAAGTTCAACTGGCCAAA; TLN-1(annealing 60°C) forward TCAAGCAGAAGTTGCACACC, reverse GGGAGCCATTAAGGATGTCA. PCR analysis was done using the ΔΔ method with 18s serving as housekeeping gene control. Statistics were done using graphpad software (PRISM version 7), paired t-test, two-tailed.

Comparison of IFNγ stimulated vs. cytomegalovirus stimulated macrophage gene expression

A total of 179 gene expression measurements were extracted from a published paper describing the fold change in expression levels of genes induced after 4 h exposure to cytomegalovirus [46]. The data was converted to a tab-delimited text file containing the official gene symbol and the reported expression level. The file was loaded into a MySQL relational database and joined to the expression data produced from the B2 and B19 cells. The data was joined on the gene symbol and a set of 54 genes were identified. The fold change for the B2 and B19 expression data was calculated by taking the log-2 (4 h expression / 0 h expression). B2 and B19 genes having expression = 0 for the initial time point were converted to 0.1 to prevent division by zero. Additionally, the fold-change reported for IL6, 280.8, was changed to 35, in order to preserve the scale of the graphs and legibility of the resulting data represented in the histograms. The fold-change in expression for the B-haplotype birds and the published data was plotted using Microsoft Excel.

Results

Differential gene expression patterns

A set of 13,618 unique genes from among all mapped sequencing reads was generated from the 4 birds across all 9 time points. Next, we analyzed the expression data to determine the number of genes expressed in each haplotype within each time point. Within the minus 6-day (t-6) time point, representing the time point after plating and adherence of monocytes and the start of differentiation into mature macrophages, 11,785 genes were expressed in the B2 birds while 12,089 were observed in the B19 birds, with 11,216 genes expressed in both. Interestingly, 4770 genes were off in both B19 and B2 haplotypes while just 569 genes exhibited expression in only the B2 chickens and 873 genes were expressed only in the B19 birds.

Similar relationships were detected in each of the remaining eight time points. The t-3 day time point, representing 3 days of differentiation in cell culture, exhibited the greatest expression of genes with a total of 11,429 expressed in both B19 and B2 birds while just 4068 genes lacked evidence of expression in both haplotypes. Also, during the t-3 day time point the greatest number of genes (1118) exhibit evidence of expression in the B2 birds while lacking evidence of expression in the B19 birds. At the t0 time point, after 6 days of differentiation and immediately before stimulation with interferon, 10,975 genes were expressed in both haplotypes while 4547 genes were not expressed in macrophages of either haplotype. Likewise, the 1 h and 2 h time points exhibited 11,349 genes and 10,789, on in both haplotypes, respectively. It is worth noting that the time point with the most genes off in both haplotypes is 16 h with 5238 genes.

Overall the data indicates that approximately 10,000 to 11,000 genes are on in both haplotypes at each time point while roughly 4000 to 5300 genes are off in both haplotypes at each time point. The number of genes on in one haplotype, while off in the other haplotype, ranges from about 400 to 1140 depending upon the haplotype and time point (Fig 1)

Fig 1

Pattern of 13,618 genes expressed across haplotypes and timepoints.

Visual representation of genes within B2 and B19 haplotypes at each of the time points. Figure includes genes expressed in common, genes expressed only in B2, genes expressed only in B19, all genes expressed in B2, all genes expressed in B19, and total non-redundant genes expressed in either B2 or B19 haplotypes.

Differences in numbers of genes expressed in B19 versus B2 haplotype birds

In order to better understand the cell biology underlying differences in macrophage differentiation and activation between B19 and B2 birds, we searched for genes exhibiting statistically significant differences between different time points within a single B-haplotype haplotype as well within the same time point between haplotypes.

When comparing the expression profiles between the B2 and B19 haplotypes, we identified 210 genes exhibiting differential expression at the t-6 day time point. These genes represent 198 genes with higher expression in the B2 birds and just 12 genes for which expression was greater in the B19. After three days, at the t-3 day time point, thousands of genes exhibited altered expression patterns between the two groups. Surprisingly, 7000 genes showed higher expression in the B19 birds while only 14 genes were expressed at higher levels in the B2 birds.

By t0 hrs, which corresponds to 6 days of monocyte differentiation into macrophages, we observed 955 genes with significant expression patterns between the haplotypes. Of these genes, 544 exhibited greater expression in the B2 haplotype while 411 exhibited higher expression in the B19 haplotype.

Cells were stimulated with IFNγ immediately following RNA collection at the t0 hr time point. At 1 h (t1) post-stimulation 665 genes show evidence of significant patterns of expression between the haplotypes where B19 birds had 109 genes expressed to higher levels while the B19 haplotype was associated with 556 genes having greater expression compared to t0. This pattern of increased expression in the B19 group is reversed by the 2 h time point.

At 2 h after IFNγ treatment, the B2 cells show a global increase in expression for 5989 genes while the B19 cells have just 18 genes on at higher levels than the B2 birds. By 4 hours after stimulation, the B2 birds still exhibit greater expression for 1029 genes while the B19 birds exhibit higher expression for 12 genes. This trend changes by 8 hours after treatment, at which time the slower responding B19 group begin showing increased expression in 797 genes while the B2 cells have greater expression for just 15 genes. By 16 hours after stimulation, only 66 genes are differentially expressed between the two haplotype groups. And, at the 24 hour mark, 406 genes show evidence of statistically significant differences in expression between them with the B2 cells exhibiting greater expression for 339 genes while the B19 cells have higher expression for 67 genes (Table 2).

Table 2

Differences in gene expression between B2 and B19.

	Total #	Number Genes	Number Genes
	Significant Genes	Higher in LEFT	Higher in RIGHT
B2 t-6 versus B19 t-6	210	198	12
B2 t-3 versus B19 t-3	7014	14	7000
B2 t0 versus B19 t0	955	544	411
B2 t1 versus B19 t1	665	109	556
B2 t2 versus B19 t2	6007	5989	18
B2 t4 versus B19 t4	1041	1029	12
B2 t8 versus B19 t8	812	15	797
B2 t16 versus B19 t16	66	28	38
B2 t24 versus B19 t24	406	339	67
B2 t-6 versus B2 t-3	6012	5998	14
B2 t-3 versus B2 t0	523	379	144
B2 t0 versus B2 t1	534	339	195
B2 t1 versus B2 t2	6104	6	6098
B2 t2 versus B2 t4	621	391	230
B2 t4 versus B2 t8	6185	6185	0
B2 t8 versus B2 t16	83	39	44
B2 t16 versus B2 t24	0	0	0
B19 t-6 versus B19 t-3	326	14	312
B19 t-3 versus B19 t0	7157	7144	13
B19 t0 versus B19 t1	67	1	66
B19 t1 versus B19 t2	180	159	21
B19 t2 versus B19 t4	1227	63	1164
B19 4 versus B19 8	70	20	50
B19 8 versus B19 16	386	362	24
B19 16 versus B19 24	24	11	13

Different temporal gene expression in B19 versus B2 haplotype birds

The B2 and B19 haplotype birds represent distinct genetic variation within the B-locus on chromosome 16. Subsequently, patterns of gene expression variation of the genes located within this region were investigated. Among the seventeen genes exhibiting statistically significant differences in expression between the B2 and B19 birds, many displayed divergent gene expression patterns prior to IFNγ stimulation. In the B2 cells, gene expression peaks on day t-6 and expression is effectively inhibited by day t-3. This is not the case in the B19 cells. Rather than reach maximum expression levels in a single day, the B19 cells don’t achieve maximum expression until day t-3 (Fig 2).

Fig 2

Distinct temporal gene expression patterns in B2 versus B19 monocytes/macrophages.

B-locus haplotypes in chickens provide a mechanism for genetically perturbing the cluster of immunologically important genes on chromosome 16 and producing phenotypic variation affecting infectious disease susceptibility and resistance. The heat map allows visualization of gene expression between the two genetically distinct haplotypes. Each row represents a gene within the B-locus (listed on the right) and each column corresponds to a particular time point when cells were collected for RNA sequencing. Black pixels indicate zero gene expression for a particular gene at a specific point in time, and dark blue corresponds to very low expression, while brighter blue indicates the next higher levels. Dark purple represents higher expression levels than blue colors, and pink represents the highest levels of gene expression. Monocytes were obtained from each haplotype of chicken and allowed to differentiate into macrophages in vitro for seven, days beginning on day minus 6 (t-6). RNA was sampled on day t-6, day t-3, and again three days later which is denoted as 0 hours (t0), when IFNγ was initially added to the cultures. On t0, RNA was sampled immediately before stimulation with IFNγ. Subsequent time points correspond to the time following interferon stimulation, in hours (1 hour, 2 hours, 4 hours, 8 hours, 16 hours and 24 hours). As visible on the heat map, there are distinct differences in gene expression between the B2 and B19 cells. The most dramatic difference occurs on day t-6. B2 cells exhibit a rapid burst of gene expression, indicated as a single column of pink on the left most edge of the heat map. In contrast, the B19 cells appear to undergo a much slower and prolonged gene expression program that was not as rapidly down regulated as in genes in the B2 cells. Additional gene expression data for a number of proteins involved in cell growth and apoptosis, is shown in the bottom half of the figure to highlight a similar pattern in gene expression and kinetics. The green border indicates the B2 haplotype expression pattern and the red border corresponds to the B19 expression pattern.

For example TRIM7, TRIM27.1, BF2, TPN, and TRIM41 exhibit strong expression on day t-6 in the B2 cells while the same genes exhibited prolonged expression over day t-6 and day t-3 in the B19 cells. Members of the TRIM (tripartite motif) family have been implicated in antiviral immune defense and several are ubiquitin ligases [47, 48]. TPN (Tapasin) is a co-factor for MHC I critical for antigen presentation to cytotoxic T-cells and chickens express the single MHCI locus termed BF-2 which is working with TPN in antigen presentation and it has been shown that there are differences in the selection of high affinity peptides in B19 vs B15 haplotypes [49] highlighting their critical role in immune competence. Additional genes within the B-locus display a similar pattern of pre-stimulatory differences in gene expression between the two different haplotypes, including genes involved in differentiation, cell growth and apoptosis such as PTPN2 (tyrosine protein phosphatase non-receptor2) and NFKB. Gene expression decreases to approximately baseline levels by time point t0 hours.

A second distinction in the gene expression patterns between B19 and B2 cells is that B2 cells exhibited a fairly robust expression at 2 and 4 hours after interferon stimulation. Unlike the B2 haplotype, the B19 haplotype appears incapable of generating such a rapid, robust and coherent gene expression profile. In contrast, the B19 cells generate a delayed, weak and uncoordinated lower level of expression that extends up to 8 hours, and in some cases even 16 hours. Overall, this global pattern of temporally dysregulated gene expression represents a re-occurring theme with the B19 monocytes and macrophages.

The divergent timing of gene expression observed in the B-locus genes is mirrored in many other genes as well, including members of the TLR signaling pathway, cellular mediators of apoptosis and cell survival, and components of cytokine signaling.

B2 and B19 display different patterns of gene expression during differentiation

The global dysregulation of gene expression among 700 genes at the t-3 day time point, as well as the expression pattern of 6000 genes exhibiting altered expression led us to explore the pattern of gene expression changes within each haplotype group over all of the time points. At the onset of the study, the B2 cells were actively expressing a diverse set of genes, however by the day t-3, most of those genes displayed reduced expression in the B2 group. Even so, the B19 haplotype cells continue to express these 7000 genes at higher levels than the B2 birds. After stimulation, B2 macrophages again show different patterns of expression compared to B19 cells in regards to timing of peak expression and coherence of expression. Four distinct patterns of divergent gene expression were identified between the B2 haplotype birds and the B19 haplotype birds (Fig 3).

Fig 3

Examples of divergent gene expression patterns observed in B2 and B19 haplotype macrophages.

Four distinct patterns were identified as representative of the types of divergent gene expression that re-occur across many genes involved in macrophage differentiation, activation and function in B2 versus B19 macrophages. 1. Day t-6: B2 high vs B19 low. This divergent pattern exhibits strong expression of genes on day -6 in the B2 birds while relatively low levels of expression are observed in the B19 birds at the same time point. Genes of interest include an adenosine receptor (P2RY12) 2. This example of divergent patterns is the single peak of day t-6 gene expression in the B2 haplotype cells compared to the prolonged multiple day expression until day t-3 in the B19 haplotype cells. Genes of interest include macrophage differentiation gene GATA, adenosine receptor A2A and macrophage podosome markers VCL and GSN. IFNγ Another interesting divergent gene expression pattern observed between the two haplotypes occurs after stimulation by IFNγ. There is a four-hour difference in peak expression timing for a large number of induced genes. In the B2 haplotype macrophages, the peak expression occurs between 2 and 4 hours, while in the B19 macrophages, the peak expression occurs between 4 and 8 hours. IFNγ Another discernable difference in post-stimulatory induction of genes between the B2 macrophages compared to the B19 macrophages is one of coherence. Specifically, there are a number of genes for which the B2 macrophages are able to rapidly turn on and reach relatively high levels of expression within 2 to 4 hours of IFNγ stimulation. In contrast, these same genes fail to exhibit a coherent peak of expression, even after 4 to 8 hours, in the B19 cells. Instead, they exhibit a dispersed “smear” of gene expression extending from approximately 1 hour after stimulation to 16 hours post-stimulation.

The first interesting divergent pattern shows strong gene expression on day t-6 in the B2 birds while relatively low levels of expression are observed in the B19 birds on the same day. This pattern is of interest because it represents a group of genes that are differentially regulated at the onset of the experimental time course. Specifically, these genes include the macrophage M1 marker PTGS2, as well as the B-locus gene cyp21. Other genes exhibiting this pattern include secreted interleukin ligands IL-1β, IL4I1, and IL6, along with genes associated with inhibition of cellular processes including IRG1 and MIP-3α. Interestingly, the adenosine receptor also displays this pattern of expression. These genes may represent initial modulators of divergent monocyte to macrophage differentiation between the B2 and B19 cells.

The second example of divergent expression patterns is the single peak of day t-6 expression in the B2 haplotype cells compared to the prolonged multiple day expression in the B19 haplotype cells. Some of these genes are macrophage differentiation mediators, like GATA2 [50], and FADD, while others are macrophage podosome (primary matrix structure) markers, including VCL and GSN. Other genes exhibiting this divergent expression pattern include chemokine receptors, like CxCR4, fatty acid transport, such as SLC25A17, and ubiquitin related factors, like DD5, which is associated with proteasomal degradation of gene products.

Additional interesting divergent gene expression patterns were observed between the two haplotypes occurring after stimulation by IFNγ (Fig 3). A notable difference in post-stimulatory induction of gene expression is a four-hour difference in peak expression timing for a large number of induced genes. In the B2 haplotype macrophages, the peak expression occurs between 2 and 4 hours, while in the B19 macrophages, the peak expression occurs between 4 and 8 hours. Some of the most noticeable genes exhibiting this divergent gene expression pattern include LITAF, IL-1β, IL12, and IFIH1, genes involved in macrophage signaling and M1 macrophage polarization [26]. Additionally, a number of genes implicated in invadosome assembly and function also exhibit this temporally displaced pattern of induction such as CD44, RAC1, and SRC.

Another discernable difference in post-stimulatory induction of genes between the B2 macrophages compared to the B19 macrophages is one of coherence (Fig 3). Specifically, there are a number of genes for which the B2 macrophages are able to rapidly turn on and reach relatively high levels of expression within 2 to 4 hours of IFNγ stimulation. In contrast, these same genes fail to exhibit a coherent peak of expression, even after 4 to 8 hours, in the B19 cells. Instead, they exhibit a dispersed “smear” of gene expression extending from approximately 1 hour after stimulation to 16 hours post-stimulation. Some of the most represented genes exhibiting this divergent pattern of expression include molecules involved in lysosome function and phagocytosis. CTTN and ACTR3, genes implicated in FcR mediated phagocytosis, along with lysosomal-associated molecules, like LAPTM5 and LAMP1, as well as the lysosomal transporter molecules ATP6AP1, ATP6V1G1 and ATP6V0C, exhibit this non-coherent pattern of expression in the B19 macrophages.

In contrast, immediately following stimulation, the B2 cells rapidly induce expression of roughly 6000 genes by the 2 h following stimulation; while, at the same time, the cells derived from the B19 birds show no signs of induction among these genes until after 4 h. It is interesting to note that while the B2 birds show a statistically significant increase in expression for 6100 genes between 1 h and 2 h, the B19 cells exhibit increased expression for just 66 genes at this time point. The largest wave of increased gene expression occurs in the B19 cells during the transition from 2 h to 4 h post stimulation, when 1164 genes increase significantly over this time period.

At the transition between 8 h and 16 h, the B2 haplotype group only exhibits differences in expression for 83 genes, with 44 having higher expression at the 16 h time point. Yet, the B19 cells show differences in 386 genes during this same period, but interestingly, 356 of these genes exhibit decreased expression during this same time interval. Taken together, these results suggest that a global disruption of temporal gene expression underlies the observed differences in differentiation, activation and nitric oxide production from macrophages derived from the two different MHC haplotypes.

RT-PCR of B2 and B19 haplotype cells following IFNγ stimulation

Gene expression was measured in separate samples of B2 and B19 cells following stimulation with IFNγ. Change in expression was assessed at 2 hours and 4 hours post stimulation. ATP6V0C exhibited the greatest induction of all genes assayed, showing an increased expression in the B2 cells at 4 hours that was 20 times the initial expression at 0-hours. Expression of ATP6VOC was dramatically less in the B19 birds. Similarly, IL18R exhibited greater than 9 times the initial expression in the B2 cells at 4 hours compared to the B19 cells which exhibited less than 2 times the initial expression at 0-hours. LITAF and TLR2 exhibited more than 7 times the expression at 4 hours in the B2 macrophages, while TLN-1, TLR-5, TLR-6 and TLR-7 exhibited greater than 4 times the initial expression in the B2 macrophages. In contrast, the B19 macrophages failed to exhibit comparable induction of these genes (Fig 4).

Fig 4

RT-PCR validation of transcripts identified as significantly expressed in RNA sequencing data.

Gene expression for ATP6V0C, LITAF, IL18R, TLN-1, TLR2, TLR3, TLR4, TLR5, TLR6, and TLR7 was assessed in B2 and B19 monocytes/macrophages following stimulation with IFNγ. Expression was measured at 0 hours, 2 hours and 4 hours. Expression for transcripts in B2 cells are shown in green and expression for transcripts are shown in red. Standard error is shown for each value. Values were considered statistically significant with p<0.05.

IFNγ stimulated vs. cytomegalovirus stimulated macrophage gene expression

In addition to the RT-PCR validation of gene expression, 54 genes, for which gene expression changes were described following cytomegalovirus stimulation were used as comparisons for the corresponding genes in the B2 and b19 haplotype birds (Fig 5). A total of 25 published genes exhibited decreases in expression following cytomegalovirus stimulation while 29 genes exhibited increased expression following stimulation. Interestingly, all but one gene (FEZ1) in the B2 cells exhibited increased expression following IFNγ stimulation. In contrast, ten genes displayed decreased expression in the B19 cells. Of the ten exhibiting fold-change < 0 in the B19 cells, 70% also exhibited decreased expression in the cytomegalovirus stimulated cells. In total, 28 genes (52%) expressed in the B2 cells matched the direction of the fold change reported in the published data while 33 genes (61%) corresponded between the B19 cells and the published data. Of the ten published genes reported as having greater than 5-fold increased expression, 90% of the B2 genes exhibited fold-change in the same direction. Overall, this data, in conjunction with the RT-PCR data, provides a comprehensive set of validation data providing evidence that the B2 and B19 gene expression data is reproducible and similar to expression patterns observed in cells stimulated towards macrophage activation pathway.

Fig 5

IFNγ stimulated vs. cytomegalovirus stimulated macrophage gene expression.

54 genes, for which gene expression changes were previously described following cytomegalovirus stimulation were used as comparisons for the corresponding genes in the B2 and B19 haplotype birds. A total of 25 published genes exhibited decreases in expression following cytomegalovirus stimulation while 29 genes exhibited increased expression following stimulation. All but one gene (FEZ1) in the B2 cells exhibited increased expression following IFNγ stimulation. In contrast, ten genes displayed decreased expression in the B19 cells. Of the ten exhibiting fold-change < 0 in the B19 cells, seven exhibited decreased expression in the cytomegalovirus stimulated cells. Twenty-eight genes (52%) expressed in the B2 cells matched the direction of the fold change reported in the published data while 33 genes (61%) corresponded between the B19 cells and the published data. Of the ten published genes reported as having greater than at least 5-fold increased expression, 90% of the B2 genes exhibited fold-change in the same direction.

Divergent non-coding RNA expression in B2 and B19 macrophages

Visualization of gene expression via heat maps facilitated the identification of distinct expression patterns between the B2 and B19 haplotypes. Because any initial differences in gene expression existing 6 days before IFNγ stimulation represent candidates responsible for the observed phenotypic differences between the two haplotypes. Genes exhibiting divergent gene expression patterns between B2 and B19 birds on day -6 were identified (Fig 6). The genes cluster into four major clades (clade1, clade2, clade3, and clade4 with a singleton labelled clade 5). Among these genes, represented in clade1 and clade2, are a number of miRNAs exhibiting strong expression in B2 cells (mir-147, mir-146b, mir-1618, mir-200a, mir-1649, and mir-1648a) compared to the B19 samples. Likewise, miRNAs contained in clade3 and clade4 exhibit greater expression in B19 cells (mir-1627, mir-222b, mir-1633, and mir-19a).

Fig 6

Identification of divergent gene expression patterns between B2 and B19 macrophages.

Visualization of divergent gene expression patterns between the B2 and B19 haplotypes. A subset of genes exhibiting divergent gene expression were identified and visualized in heat following hierarchical clustering of the genes (rows), but not the time points (columns). The genes cluster into four major clades (clade1, clade2, clade3, and clade4) with a singleton gene (labelled clade 5). Among these genes, represented in clade1 and clade2, are a number of miRNAs exhibiting strong expression in B2 cells (mir-147, mir-146b, mir-1618, mir-200a, mir-1649, and mir-1648a) compared to the B19 samples. Likewise, miRNAs contained in clade3 and clade4 exhibit greater expression in B19 cells (mir-1627, mir-222b, mir-1633, and mir-19a). Additionally, a number of small nucleolar RNAs (snoRNAs) exhibit similarly dichotomous gene expression patterns (clade4) such that SNORd24, snoZ40, SNORD74, SNORA17, and SNORD12 exhibit substantially higher levels of expression in B19 cells on day -6 compared to B2 cells while B2 cells express such as snoU2_19 (clade2).

A number of small nucleolar RNAs (snoRNAs) exhibit similarly dichotomous gene expression patterns (clade4). For example, SNORd24, snoZ40, SNORD74, SNORA17, and SNORD12 exhibit substantially higher levels of expression in B19 cells on day -6 compared to B2 cells (Fig 6). However, B2 cells also express snoRNAs exhibiting divergent expression patterns between the two haplotypes, such as snoU2_19 (clade2).

In addition to non-coding RNAs, divergent expression patterns are also observed with protein-coding RNAs (Fig 6). For example, clade1 contains IL6, IL18, IL-1β, CCL1, PTPN2 and MMP10, which exhibit higher initial expression on day -6 in the B2 birds. In contrast, the protein-coding genes LAMP2, UBXN7, UBE4B, PK3CA, UBE2W, and CX3CR1, in clade3, exhibit higher initial expression patterns in B19 birds. Clade5 contains the single gene IFNγ, which exhibits relatively low expression early in both B2 and B19 cells, but following stimulation rises to a higher level at 8 hours in the B19 birds.

Considering the diverse expression patterns discovered and the results indicating the involvement of non-coding RNAs, further results including divergent non-coding RNA expression will be described in more detail in further publications.

Gene enrichment analysis

Enrichment analysis of genes exhibiting statistically significant differences in expression between time points and/or haplotypes (Table 3, S2, S3 and S4 Tables) was performed using gene ontology and both KEGG and reactome pathways. The results of the gene enrichment provided a high-resolution perspective of the functional role of mRNA sequenced within the B2 cells across the experimental time points.

Table 3

Gene enrichment analysis—Highlights from gene ontology, KEGG pathways, reactome pathways.

Sample Comparison	Gene Set	Term	Description	Gene Count	P-Value
B2 t-6 vs. B2 t-3	Down	GO:0006360	transcription from RNA polymerase I promoter	12	8.08E-05
B2 t-6 vs. B2 t-3	Down	GO:0000398	mRNA splicing, via spliceosome	33	8.85E-05
B2 t-6 vs. B2 t-3	Down	GO:0043966	histone H3 acetylation	20	4.18E-04
B2 t-6 vs. B2 t-3	Down	GO:0008333	endosome to lysosome transport	17	0.001611612
B2 t-6 vs. B2 t-3	Down	GO:0008033	tRNA processing	12	0.002902509
B2 t-6 vs. B2 t-3	Down	GO:0031338	regulation of vesicle fusion	15	0.005935057
B2 t-6 vs. B2 t-3	Down	GO:0006378	mRNA polyadenylation	11	0.006175743
B2 t-6 vs. B2 t-3	Down	GO:0000387	spliceosomal snRNP assembly	13	0.006325454
B2 t-6 vs. B2 t-3	Down	GO:0006397	mRNA processing	27	0.010374791
B2 t-6 vs. B2 t-3	Down	GO:0006886	intracellular protein transport	61	0.012613014
B2 t-6 vs. B2 t-3	Down	GO:0043123	positive regulation of I-kappaB kinase/NF-kappaB signaling	42	0.013554129
B2 t-6 vs. B2 t-3	Down	GO:0016050	vesicle organization	11	0.023348647
B2 t-6 vs. B2 t-3	Down	GO:0030968	endoplasmic reticulum unfolded protein response	16	0.023521606
B2 t-6 vs. B2 t-3	Down	GO:0045022	early endosome to late endosome transport	10	0.024838372
B2 t-6 vs. B2 t-3	Down	GO:0045292	mRNA cis splicing, via spliceosome	7	0.025727402
B2 t-6 vs. B2 t-3	Down	GO:0016226	iron-sulfur cluster assembly	9	0.047036045
B2 t-6 vs. B2 t-3	Down	GO:0007032	endosome organization	13	0.048431064
B2 t-6 vs. B2 t-3	Down	GO:0032088	negative regulation of NF-kappaB transcription factor activity	16	0.048807155
B2 t-6 vs. B2 t-3	Down	gga03020	RNA polymerase	17	8.84E-05
B2 t-6 vs. B2 t-3	Down	gga03018	RNA degradation	38	6.00E-04
B2 t-6 vs. B2 t-3	Down	gga01100	Metabolic pathways	432	0.001331836
B2 t-6 vs. B2 t-3	Down	gga03008	Ribosome biogenesis in eukaryotes	37	0.002017081
B2 t-6 vs. B2 t-3	Down	gga03040	Spliceosome	54	0.002180456
B2 t-6 vs. B2 t-3	Down	gga00190	Oxidative phosphorylation	57	0.004584223
B2 t-6 vs. B2 t-3	Down	gga03022	Basal transcription factors	21	0.021845672
B2 t-6 vs. B2 t-3	Down	gga03060	Protein export	14	0.025867798
B2 t-6 vs. B2 t-3	Down	gga03015	mRNA surveillance pathway	34	0.034772067
B2 t-6 vs B2 t-3	Down	R-GGA-5419276	Mitochondrial translation termination	41	4.06E-10
B2 t-6 vs B2 t-3	Down	R-GGA-5389840	Mitochondrial translation elongation	39	5.57E-09
B2 t-6 vs B2 t-3	Down	R-GGA-73779	RNA Polymerase II Transcription Pre-Initiation And Promoter Opening	24	5.04E-06
B2 t-6 vs B2 t-3	Down	R-GGA-75953	RNA Polymerase II Transcription Initiation	24	5.04E-06
B2 t-6 vs B2 t-3	Down	R-GGA-76042	RNA Polymerase II Transcription Initiation And Promoter Clearance	24	5.04E-06
B2 t-6 vs B2 t-3	Down	R-GGA-674695	RNA Polymerase II Pre-transcription Events	34	5.58E-06
B2 t-6 vs B2 t-3	Down	R-GGA-72086	mRNA Capping	19	1.15E-05
B2 t-6 vs B2 t-3	Down	R-GGA-75955	RNA Polymerase II Transcription Elongation	24	1.38E-04
B2 t-6 vs B2 t-3	Down	R-GGA-72165	mRNA Splicing—Minor Pathway	26	2.60E-04
B2 t-6 vs B2 t-3	Down	R-GGA-72163	mRNA Splicing—Major Pathway	51	6.17E-04
B2 t-6 vs B2 t-3	Down	R-GGA-983168	Antigen processing: Ubiquitination & Proteasome degradation	32	0.001262633
B2 t-6 vs B2 t-3	Down	R-GGA-1834949	Cytosolic sensors of pathogen-associated DNA	12	0.001290237
B2 t-6 vs B2 t-3	Down	R-GGA-611105	Respiratory electron transport	24	0.002379624
B2 t-6 vs B2 t-3	Down	R-GGA-1855183	Synthesis of IP2, IP, and Ins in the cytosol	7	0.009896133
B2 t-6 vs B2 t-3	Down	R-GGA-180292	GAB1 signalosome	8	0.028654133
B2 t-6 vs B2 t-3	Down	R-GGA-189451	Heme biosynthesis	7	0.027634144
B2 t-3 vs. B2 t0	Down	GO:0007059	chromosome segregation	8	1.87E-05
B2 t-3 vs. B2 t0	Down	GO:0007067	mitotic nuclear division	7	0.002285326
B2 t-3 vs. B2 t0	Down	GO:0007018	microtubule-based movement	6	0.003424586
B2 t-3 vs. B2 t0	Down	GO:0000281	mitotic cytokinesis	4	0.004035208
B2 t-3 vs. B2 t0	Down	GO:0008152	metabolic process	6	0.01138645
B2 t-3 vs. B2 t0	Down	GO:0006281	DNA repair	7	0.021563304
B2 t-3 vs. B2 t0	Down	GO:0045671	negative regulation of osteoclast differentiation	3	0.02255698
B2 t-3 vs. B2 t0	Down	GO:0051301	cell division	6	0.030134621
B2 t-3 vs. B2 t0	Down	gga04110	Cell cycle	14	5.77E-06
B2 t-3 vs. B2 t0	Down	gga03030	DNA replication	7	1.16E-04
B2 t-3 vs. B2 t0	Down	gga03430	Mismatch repair	5	0.00172114
B2 t-3 vs. B2 t0	Down	gga00240	Pyrimidine metabolism	9	0.002460978
B2 t-3 vs. B2 t0	Down	gga04630	Jak-STAT signaling pathway	8	0.034070609
B2 t-3 vs. B2 t0	Down	R-GGA-5663220	RHO GTPases Activate Formins	14	1.85E-07
B2 t-3 vs. B2 t0	Down	R-GGA-2500257	Resolution of Sister Chromatid Cohesion	14	2.19E-07
B2 t-3 vs. B2 t0	Down	R-GGA-2467813	Separation of Sister Chromatids	14	1.41E-06
B2 t-3 vs. B2 t0	Down	R-GGA-69205	G1/S-Specific Transcription	5	2.53E-04
B2 t-3 vs. B2 t0	Down	R-GGA-113510	E2F mediated regulation of DNA replication	5	2.53E-04
B2 t-3 vs. B2 t0	Down	R-GGA-983189	Kinesins	4	0.001487779
B2 t-3 vs. B2 t0	Down	R-GGA-156582	Acetylation	3	0.002846986
B2 t-3 vs. B2 t0	Down	R-GGA-5358565	Mismatch repair (MMR) directed by MSH2:MSH6 (MutSalpha)	4	0.003044778
B2 t-3 vs. B2 t0	Down	R-GGA-5651801	PCNA-Dependent Long Patch Base Excision Repair	4	0.00409159
B2 t-3 vs. B2 t0	Down	R-GGA-512988	Interleukin-3, 5 and GM-CSF signaling	4	0.006774634
B2 t-3 vs. B2 t0	Down	R-GGA-912526	Interleukin receptor SHC signaling	3	0.030200217
B2 t-3 vs. B2 t0	Up	GO:0071353	cellular response to interleukin-4	4	1.26E-04
B2 t-3 vs. B2 t0	Up	GO:0006564	L-serine biosynthetic process	2	0.025467339
B2 t-3 vs. B2 t0	Up	GO:0006166	purine ribonucleoside salvage	2	0.033813353
B2 t-3 vs. B2 t0	Up	GO:0006366	transcription from RNA polymerase II promoter	4	0.047682647
B2 t-3 vs. B2 t0	Up	gga04141	Protein processing in endoplasmic reticulum	7	0.001678292
B2 t-3 vs. B2 t0	Up	gga01230	Biosynthesis of amino acids	4	0.015296859
B2 t-3 vs. B2 t0	Up	gga01100	Metabolic pathways	16	0.036140022
B2 t-3 vs. B2 t0	Up	gga00260	Glycine, serine and threonine metabolism	3	0.03864845
B2 t-3 vs. B2 t0	Up	R-GGA-977347	Serine biosynthesis	2	0.023430989
B2 t-3 vs. B2 t0	Up	R-GGA-433692	Proton-coupled monocarboxylate transport	2	0.046329225
B2 t0 vs B2 t1	Down	GO:0006412	translation	17	2.85E-08
B2 t0 vs B2 t1	Down	GO:0042149	cellular response to glucose starvation	6	3.11E-05
B2 t0 vs B2 t1	Down	GO:0006457	protein folding	9	4.18E-04
B2 t0 vs B2 t1	Down	GO:0030968	endoplasmic reticulum unfolded protein response	4	0.023705124
B2 t0 vs B2 t1	Down	GO:0030970	retrograde protein transport, ER to cytosol	3	0.02551737
B2 t0 vs. B2 t1	Down	gga03010	Ribosome	21	7.21E-11
B2 t0 vs. B2 t1	Down	gga04141	Protein processing in endoplasmic reticulum	20	2.24E-08
B2 t0 vs. B2 t1	Down	gga00970	Aminoacyl-tRNA biosynthesis	7	0.001164589
B2 t0 vs. B2 t1	Down	gga00330	Arginine and proline metabolism	6	0.006099731
B2 t0 vs. B2 t1	Down	R-GGA-1799339	SRP-dependent cotranslational protein targeting to membrane	15	4.23E-11
B2 t0 vs. B2 t1	Down	R-GGA-72706	GTP hydrolysis and joining of the 60S ribosomal subunit	14	2.46E-10
B2 t0 vs. B2 t1	Down	R-GGA-975956	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	14	6.02E-10
B2 t0 vs. B2 t1	Down	R-GGA-975957	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	14	5.10E-09
B2 t0 vs. B2 t1	Down	R-GGA-72695	Formation of the ternary complex, and subsequently, the 43S complex	9	7.35E-07
B2 t0 vs B2 t1	Up	GO:0006954	inflammatory response	11	6.00E-06
B2 t0 vs B2 t1	Up	GO:0051607	defense response to virus	6	6.95E-04
B2 t0 vs B2 t1	Up	GO:0002224	toll-like receptor signaling pathway	4	0.002646723
B2 t0 vs B2 t1	Up	GO:0060326	cell chemotaxis	4	0.006725072
B2 t0 vs B2 t1	Up	GO:0007596	blood coagulation	4	0.012276628
B2 t0 vs B2 t1	Up	GO:0002755	MyD88-dependent toll-like receptor signaling pathway	3	0.012708351
B2 t0 vs B2 t1	Up	GO:0071222	cellular response to lipopolysaccharide	4	0.013230968
B2 t0 vs B2 t1	Up	GO:0007052	mitotic spindle organization	3	0.014697099
B2 t0 vs B2 t1	Up	GO:0006955	immune response	6	0.015620183
B2 t0 vs B2 t1	Up	GO:0002548	monocyte chemotaxis	3	0.019044739
B2 t0 vs B2 t1	Up	GO:0009263	deoxyribonucleotide biosynthetic process	2	0.03982406
B2 t0 vs. B2 t1	Up	gga04620	Toll-like receptor signaling pathway	7	0.001706356
B2 t0 vs. B2 t1	Up	gga04630	Jak-STAT signaling pathway	7	0.008611879
B2 t0 vs. B2 t1	Up	gga04514	Cell adhesion molecules (CAMs)	6	0.026832539
B2 t0 vs. B2 t1	Up	gga04068	FoxO signaling pathway	6	0.038420138
B2 t0 vs. B2 t1	Up	gga05164	Influenza A	6	0.04852434
B2 t0 vs. B2 t1	Up	gga00240	Pyrimidine metabolism	5	0.049377943
B2 t0 vs. B2 t1	Down	R-GGA-1799339	SRP-dependent cotranslational protein targeting to membrane	15	4.23E-11
B2 t0 vs. B2 t1	Down	R-GGA-72706	GTP hydrolysis and joining of the 60S ribosomal subunit	14	2.46E-10
B2 t0 vs. B2 t1	Down	R-GGA-975956	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	14	6.02E-10
B2 t0 vs. B2 t1	Down	R-GGA-975957	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	14	5.10E-09
B2 t0 vs. B2 t1	Down	R-GGA-72695	Formation of the ternary complex, and subsequently, the 43S complex	9	7.35E-07
B2 t0 vs. B2 t1	Down	R-GGA-72702	Ribosomal scanning and start codon recognition	9	1.44E-06
B2 t0 vs. B2 t1	Down	R-GGA-156590	Glutathione conjugation	3	0.020307126
B2 t0 vs. B2 t1	Down	R-GGA-5673000	RAF activation	3	0.020307126
B2 t0 vs. B2 t1	Down	R-GGA-70614	Amino acid synthesis and interconversion (transamination)	3	0.033683435
B2 t1 vs. B2 t2	Up	GO:0006511	ubiquitin-dependent protein catabolic process	54	5.01E-07
B2 t1 vs. B2 t2	Up	GO:0006886	intracellular protein transport	78	7.15E-06
B2 t1 vs. B2 t2	Up	GO:0006888	ER to Golgi vesicle-mediated transport	31	1.65E-05
B2 t1 vs. B2 t2	Up	GO:0000398	mRNA splicing, via spliceosome	35	8.42E-05
B2 t1 vs. B2 t2	Up	GO:0050821	protein stabilization	42	9.33E-05
B2 t1 vs. B2 t2	Up	GO:0045454	cell redox homeostasis	33	1.41E-04
B2 t1 vs. B2 t2	Up	GO:0007030	Golgi organization	33	6.86E-04
B2 t1 vs. B2 t2	Up	GO:0043001	Golgi to plasma membrane protein transport	13	0.001353221
B2 t1 vs. B2 t2	Up	GO:0008333	endosome to lysosome transport	18	0.001409269
B2 t1 vs. B2 t2	Up	GO:0006338	chromatin remodeling	23	0.005243778
B2 t1 vs. B2 t2	Up	GO:0019827	stem cell population maintenance	21	0.005362005
B2 t1 vs. B2 t2	Up	GO:0000920	cell separation after cytokinesis	11	0.005858508
B2 t1 vs. B2 t2	Up	GO:0071353	cellular response to interleukin-4	10	0.011871699
B2 t1 vs. B2 t2	Up	GO:0016050	vesicle organization	12	0.014441064
B2 t1 vs. B2 t2	Up	GO:0043966	histone H3 acetylation	18	0.014875667
B2 t1 vs. B2 t2	Up	GO:0000381	regulation of alternative mRNA splicing, via spliceosome	14	0.015404037
B2 t1 vs. B2 t2	Up	GO:0031338	regulation of vesicle fusion	15	0.015497515
B2 t1 vs. B2 t2	Up	GO:0034067	protein localization to Golgi apparatus	7	0.015640559
B2 t1 vs. B2 t2	Up	GO:0006606	protein import into nucleus	20	0.021272498
B2 t1 vs. B2 t2	Up	GO:0030970	retrograde protein transport, ER to cytosol	9	0.023542072
B2 t1 vs. B2 t2	Up	GO:0031398	positive regulation of protein ubiquitination	17	0.024088893
B2 t1 vs. B2 t2	Up	gga04141	Protein processing in endoplasmic reticulum	93	1.04E-07
B2 t1 vs. B2 t2	Up	gga03010	Ribosome	76	9.85E-07
B2 t1 vs. B2 t2	Up	gga04120	Ubiquitin mediated proteolysis	77	4.09E-06
B2 t1 vs. B2 t2	Up	gga03040	Spliceosome	66	1.88E-05
B2 t1 vs. B2 t2	Up	gga03018	RNA degradation	44	5.23E-05
B2 t1 vs. B2 t2	Up	gga03020	RNA polymerase	18	8.51E-05
B2 t1 vs. B2 t2	Up	gga00190	Oxidative phosphorylation	68	2.07E-04
B2 t1 vs. B2 t2	Up	gga03013	RNA transport	75	3.69E-04
B2 t1 vs. B2 t2	Up	gga03420	Nucleotide excision repair	26	6.71E-04
B2 t1 vs. B2 t2	Up	gga03060	Protein export	18	7.34E-04
B2 t1 vs. B2 t2	Up	gga00240	Pyrimidine metabolism	53	9.51E-04
B2 t1 vs. B2 t2	Up	gga00510	N-Glycan biosynthesis	30	0.002686839
B2 t1 vs. B2 t2	Up	gga04110	Cell cycle	60	0.006890719
B2 t1 vs. B2 t2	Up	gga04142	Lysosome	60	0.006890719
B2 t1 vs. B2 t2	Up	gga00071	Fatty acid degradation	21	0.021537199
B2 t1 vs. B2 t2	Up	gga03022	Basal transcription factors	23	0.021579106
B2 t1 vs. B2 t2	Up	gga03015	mRNA surveillance pathway	38	0.029025863
B2 t1 vs. B2 t2	Up	gga04144	Endocytosis	113	0.032927988
B2 t1 vs. B2 t2	Up	gga04150:	mTOR signaling pathway	28	0.039075053
B2 t1 vs. B2 t2	Up	R-GGA-5419276	Mitochondrial translation termination	40	3.56E-07
B2 t1 vs. B2 t2	Up	R-GGA-5389840	Mitochondrial translation elongation	39	6.66E-07
B2 t1 vs. B2 t2	Up	R-GGA-2467813	Separation of Sister Chromatids	57	2.33E-06
B2 t1 vs. B2 t2	Up	R-GGA-2500257:	Resolution of Sister Chromatid Cohesion	50	4.11E-06
B2 t1 vs. B2 t2	Up	R-GGA-72165	mRNA Splicing—Minor Pathway	28	4.11E-04
B2 t1 vs. B2 t2	Up	R-GGA-72086	mRNA Capping	18	6.97E-04
B2 t1 vs. B2 t2	Up	R-GGA-76042	RNA Polymerase II Transcription Initiation And Promoter Clearance	24	9.36E-05
B2 t1 vs. B2 t2	Up	R-GGA-75953	RNA Polymerase II Transcription Initiation	24	9.36E-05
B2 t1 vs. B2 t2	Up	R-GGA-73779	RNA Polymerase II Transcription Pre-Initiation And Promoter Opening	24	9.36E-05
B2 t1 vs. B2 t2	Up	R-GGA-1834949	Cytosolic sensors of pathogen-associated DNA	13	8.95E-04
B2 t1 vs. B2 t2	Up	R-GGA-166208	mTORC1-mediated signaling	10	0.007874037
B2 t1 vs. B2 t2	Up	R-GGA-5674135	MAP2K and MAPK activation	11	0.010508603
B2 t1 vs. B2 t2	Up	R-GGA-202424	Downstream TCR signaling	13	0.02759075
B2 t1 vs. B2 t2	Up	R-GGA-2871796	FCERI mediated MAPK activation	14	0.029000225
B2 t1 vs. B2 t2	Up	R-GGA-174084	Autodegradation of Cdh1 by Cdh1:APC/C	24	0.029609758
B2 t1 vs. B2 t2	Up	R-GGA-2730905	Role of LAT2/NTAL/LAB on calcium mobilization	8	0.031070119
B2 t1 vs. B2 t2	Up	R-GGA-5607764	CLEC7A (Dectin-1) signaling	10	0.041042518
B2 t2 vs. B2 t4	Down	GO:0090307	mitotic spindle assembly	5	6.17E-04
B2 t2 vs. B2 t4	Down	GO:0000070	mitotic sister chromatid segregation	4	0.001514614
B2 t2 vs. B2 t4	Down	GO:0007059	chromosome segregation	5	0.004071808
B2 t2 vs. B2 t4	Down	GO:0007094	mitotic spindle assembly checkpoint	3	0.008710879
B2 t2 vs. B2 t4	Down	GO:0046849	bone remodeling	3	0.011065539
B2 t2 vs. B2 t4	Down	GO:0006464	cellular protein modification process	3	0.013666474
B2 t2 vs. B2 t4	Down	GO:0035556	intracellular signal transduction	10	0.014298815
B2 t2 vs. B2 t4	Down	gga04110	Cell cycle	11	3.28E-05
B2 t2 vs. B2 t4	Down	R-GGA-2467813	Separation of Sister Chromatids	11	2.41E-05
B2 t2 vs. B2 t4	Down	R-GGA-5663220	RHO GTPases Activate Formins	10	3.80E-05
B2 t2 vs. B2 t4	Down	R-GGA-2500257	Resolution of Sister Chromatid Cohesion	10	4.26E-05
B2 t2 vs. B2 t4	Down	R-GGA-5620912	Anchoring of the basal body to the plasma membrane	6	0.016443642
B2 t2 vs. B2 t4	Down	R-GGA-5693568	Resolution of D-loop Structures through Holliday Junction Intermediates	4	0.016488095
B2 t2 vs. B2 t4	Down	R-GGA-5620922	BBSome-mediated cargo-targeting to cilium	3	0.022369863
B2 t2 vs. B2 t4	Down	R-GGA-69205	G1/S-Specific Transcription	3	0.026920963
B2 t2 vs. B2 t4	Down	R-GGA-113510	E2F mediated regulation of DNA replication	3	0.026920963
B2 t2 vs. B2 t4	Down	R-GGA-176187	Activation of ATR in response to replication stress	4	0.027999251
B2 t2 vs. B2 t4	Down	R-GGA-606279	Deposition of new CENPA-containing nucleosomes at the centromere	4	0.03650445
B2 t2 vs. B2 t4	Down	R-GGA-2565942	Regulation of PLK1 Activity at G2/M Transition	5	0.038794821
B2 t2 vs. B2 t4	Up	GO:0006955	immune response	12	1.35E-07
B2 t2 vs. B2 t4	Up	GO:0006954	inflammatory response	11	4.10E-06
B2 t2 vs. B2 t4	Up	GO:0032496	response to lipopolysaccharide	8	4.67E-06
B2 t2 vs. B2 t4	Up	GO:0042832	defense response to protozoan	3	0.004389284
B2 t2 vs. B2 t4	Up	GO:0051607	defense response to virus	5	0.004729744
B2 t2 vs. B2 t4	Up	GO:0097190	apoptotic signaling pathway	4	0.004879708
B2 t2 vs. B2 t4	Up	GO:0042981	regulation of apoptotic process	6	0.007547916
B2 t2 vs. B2 t4	Up	GO:0032735	positive regulation of interleukin-12 production	3	0.008406525
B2 t2 vs. B2 t4	Up	GO:0045071	negative regulation of viral genome replication	3	0.008406525
B2 t2 vs. B2 t4	Up	GO:0042742	defense response to bacterium	4	0.008606681
B2 t2 vs. B2 t4	Up	GO:0043410	positive regulation of MAPK cascade	4	0.009352126
B2 t2 vs. B2 t4	Up	GO:0042127	regulation of cell proliferation	6	0.011369564
B2 t2 vs. B2 t4	Up	GO:0050717	positive regulation of interleukin-1 alpha secretion	2	0.025630045
B2 t2 vs. B2 t4	Up	GO:0048873	homeostasis of number of cells within a tissue	3	0.026931634
B2 t2 vs. B2 t4	Up	gga04060	Cytokine-cytokine receptor interaction	13	5.07E-07
B2 t2 vs. B2 t4	Up	gga05164	Influenza A	7	0.004668633
B2 t4 vs. B2 t8	Down	GO:0042787	protein ubiquitination involved in ubiquitin-dependent protein catabolic process	63	1.69E-07
B2 t4 vs. B2 t8	Down	GO:0006511	ubiquitin-dependent protein catabolic process	55	9.75E-07
B2 t4 vs. B2 t8	Down	GO:0030433	ER-associated ubiquitin-dependent protein catabolic process	30	3.12E-06
B2 t4 vs. B2 t8	Down	GO:0006886	intracellular protein transport	81	5.55E-06
B2 t4 vs. B2 t8	Down	GO:0045454	cell redox homeostasis	35	4.20E-05
B2 t4 vs. B2 t8	Down	GO:0006888	ER to Golgi vesicle-mediated transport	31	5.01E-05
B2 t4 vs. B2 t8	Down	GO:0007030	Golgi organization	36	8.73E-05
B2 t4 vs. B2 t8	Down	GO:0050821	protein stabilization	43	1.27E-04
B2 t4 vs. B2 t8	Down	GO:0000398	mRNA splicing, via spliceosome	35	2.60E-04
B2 t4 vs. B2 t8	Down	GO:0006457	protein folding	48	3.44E-04
B2 t4 vs. B2 t8	Down	GO:0006360	transcription from RNA polymerase I promoter	12	3.71E-04
B2 t4 vs. B2 t8	Down	GO:0000209	protein polyubiquitination	37	4.73E-04
B2 t4 vs. B2 t8	Down	GO:0015031	protein transport	43	5.18E-04
B2 t4 vs. B2 t8	Down	GO:0008333	endosome to lysosome transport	19	6.62E-04
B2 t4 vs. B2 t8	Down	GO:0043161	proteasome-mediated ubiquitin-dependent protein catabolic process	48	7.64E-04
B2 t4 vs. B2 t8	Down	GO:0019827	stem cell population maintenance	21	0.010108895
B2 t4 vs. B2 t8	Down	GO:0007049	cell cycle	26	0.010146854
B2 t4 vs. B2 t8	Down	GO:0031929	TOR signaling	10	0.017015043
B2 t4 vs. B2 t8	Down	gga04141	Protein processing in endoplasmic reticulum	98	8.59E-09
B2 t4 vs. B2 t8	Down	gga03010	Ribosome	78	9.30E-07
B2 t4 vs. B2 t8	Down	gga04142	Lysosome	73	1.78E-06
B2 t4 vs. B2 t8	Down	gga03040	Spliceosome	69	6.04E-06
B2 t4 vs. B2 t8	Down	gga00190	Oxidative phosphorylation	72	4.26E-05
B2 t4 vs. B2 t8	Down	gga03013	RNA transport	77	4.14E-04
B2 t4 vs. B2 t8	Down	gga04120	Ubiquitin mediated proteolysis	72	7.75E-04
B2 t4 vs. B2 t8	Down	gga03020	RNA polymerase	17	9.00E-04
B2 t4 vs. B2 t8	Down	gga03050	Proteasome	26	0.001325244
B2 t4 vs. B2 t8	Down	gga00510	N-Glycan biosynthesis	31	0.002147473
B2 t4 vs. B2 t8	Down	gga03018	RNA degradation	41	0.002187297
B2 t4 vs. B2 t8	Down	gga03015	mRNA surveillance pathway	43	0.002311437
B2 t4 vs. B2 t8	Down	gga00071	Fatty acid degradation	23	0.005623955
B2 t4 vs. B2 t8	Down	gga04150	mTOR signaling pathway	30	0.017796114
B2 t4 vs. B2 t8	Down	gga04144	Endocytosis	118	0.023369828
B2 t4 vs. B2 t8	Down	gga00562	Inositol phosphate metabolism	39	0.023965717
B2 t4 vs. B2 t8	Down	gga00240	Pyrimidine metabolism	49	0.026875801
B2 t4 vs. B2 t8	Down	gga01100	Metabolic pathways	489	0.030757983
B2 t4 vs. B2 t8	Down	gga03022	Basal transcription factors	23	0.03456006
B2 t4 vs. B2 t8	Down	gga04621	NOD-like receptor signaling pathway	25	0.03499728
B2 t4 vs. B2 t8	Down	gga04068	FoxO signaling pathway	63	0.03789817
B2 t4 vs. B2 t8	Down	R-GGA-73762	RNA Polymerase I Transcription Initiation	29	4.07E-08
B2 t4 vs. B2 t8	Down	R-GGA-72163	mRNA Splicing—Major Pathway	68	6.37E-08
B2 t4 vs. B2 t8	Down	R-GGA-5419276	Mitochondrial translation termination	41	2.05E-07
B2 t4 vs. B2 t8	Down	R-GGA-73772	RNA Polymerase I Promoter Escape	22	1.56E-06
B2 t4 vs. B2 t8	Down	R-GGA-5389840	Mitochondrial translation elongation	39	1.73E-06
B2 t4 vs. B2 t8	Down	R-GGA-674695	RNA Polymerase II Pre-transcription Events	38	3.11E-06
B2 t4 vs. B2 t8	Down	R-GGA-1799339	SRP-dependent cotranslational protein targeting to membrane	45	6.11E-06
B2 t4 vs. B2 t8	Down	R-GGA-73863	RNA Polymerase I Transcription Termination	21	1.53E-05
B2 t4 vs. B2 t8	Down	R-GGA-6781823	Formation of TC-NER Pre-Incision Complex	32	1.93E-05
B2 t4 vs. B2 t8	Down	R-GGA-975957	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	50	2.13E-05
B2 t4 vs. B2 t8	Down	R-GGA-73779	RNA Polymerase II Transcription Pre-Initiation And Promoter Opening	25	3.47E-05
B2 t4 vs. B2 t8	Down	R-GGA-75953	RNA Polymerase II Transcription Initiation	25	3.47E-05
B2 t4 vs. B2 t8	Down	R-GGA-76042	RNA Polymerase II Transcription Initiation And Promoter Clearance	25	3.47E-05
B2 t4 vs. B2 t8	Down	R-GGA-72706	GTP hydrolysis and joining of the 60S ribosomal subunit	41	4.75E-05
B2 t4 vs. B2 t8	Down	R-GGA-975956	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	43	5.97E-05
B2 t4 vs. B2 t8	Down	R-GGA-75955	RNA Polymerase II Transcription Elongation	27	6.68E-05
B2 t4 vs. B2 t8	Down	R-GGA-5696395	Formation of Incision Complex in GG-NER	25	9.14E-05
B2 t4 vs. B2 t8	Down	R-GGA-6782210	Gap-filling DNA repair synthesis and ligation in TC-NER	32	9.24E-05
B2 t4 vs. B2 t8	Down	R-GGA-6782135	Dual incision in TC-NER	33	1.10E-04
B2 t4 vs. B2 t8	Down	R-GGA-72086	mRNA Capping	19	2.11E-04
B2 t4 vs. B2 t8	Down	R-GGA-72165	mRNA Splicing—Minor Pathway	28	7.54E-04
B2 t4 vs. B2 t8	Down	R-GGA-204005	COPII (Coat Protein 2) Mediated Vesicle Transport	30	0.006870202
B2 t4 vs. B2 t8	Down	R-GGA-1834949	Cytosolic sensors of pathogen-associated DNA	12	0.007282044
B2 t4 vs. B2 t8	Down	R-GGA-983168	Antigen processing: Ubiquitination & Proteasome degradation	34	0.007438297
B2 t4 vs. B2 t8	Down	R-GGA-202424	Downstream TCR signaling	14	0.011275517
B2 t4 vs. B2 t8	Down	R-GGA-917729	Endosomal Sorting Complex Required For Transport (ESCRT)	14	0.011275517
B2 t4 vs. B2 t8	Down	R-GGA-2730905	Role of LAT2/NTAL/LAB on calcium mobilization	8	0.037401415
B2 t4 vs. B2 t8	Down	R-GGA-2871796	FCERI mediated MAPK activation	14	0.038198199
B2 t4 vs. B2 t8	Down	R-GGA-109688	Cleavage of Growing Transcript in the Termination Region	23	0.043691582
B2 t4 vs. B2 t8	Down	R-GGA-166208	mTORC1-mediated signaling	9	0.044990961
B2 t4 vs. B2 t8	Down	R-GGA-1445148	Translocation of GLUT4 to the plasma membrane	9	0.044990961
B2 t8 vs. B2 t16	Down	GO:0010634	positive regulation of epithelial cell migration	4	7.45E-05
B2 t8 vs. B2 t16	Down	GO:0045747	positive regulation of Notch signaling pathway	4	2.45E-04
B2 t8 vs. B2 t16	Down	GO:0006955	immune response	6	3.59E-04
B2 t8 vs. B2 t16	Down	GO:0006954	inflammatory response	6	6.50E-04
B2 t8 vs. B2 t16	Down	GO:0032735	positive regulation of interleukin-12 production	3	0.001678191
B2 t8 vs. B2 t16	Down	GO:0040008	regulation of growth	3	0.002746603
B2 t8 vs. B2 t16	Down	GO:0051607	defense response to virus	4	0.00329606
B2 t8 vs. B2 t16	Down	GO:2000379	positive regulation of reactive oxygen species metabolic process	3	0.004551512
B2 t8 vs. B2 t16	Down	GO:0018107	peptidyl-threonine phosphorylation	3	0.008034447
B2 t8 vs. B2 t16	Down	GO:0090002	establishment of protein localization to plasma membrane	3	0.010846074
B2 t8 vs. B2 t16	Down	GO:0018401	peptidyl-proline hydroxylation to 4-hydroxy-L-proline	2	0.016916932
B2 t8 vs. B2 t16	Down	GO:0007050	cell cycle arrest	3	0.017563794
B2 t8 vs. B2 t16	Down	GO:2000107	negative regulation of leukocyte apoptotic process	2	0.022493312
B2 t8 vs. B2 t16	Down	GO:0007179	transforming growth factor beta receptor signaling pathway	3	0.026723558
B2 t8 vs. B2 t16	Down	GO:0070102	interleukin-6-mediated signaling pathway	2	0.028038678
B2 t8 vs. B2 t16	Down	GO:2000505	regulation of energy homeostasis	2	0.033553199
B2 t8 vs. B2 t16	Down	GO:0009612	response to mechanical stimulus	2	0.033553199
B2 t8 vs. B2 t16	Down	GO:0032496	response to lipopolysaccharide	3	0.036148438
B2 t8 vs. B2 t16	Down	gga04060	Cytokine-cytokine receptor interaction	9	2.22E-05
B2 t8 vs. B2 t16	Down	gga04630	Jak-STAT signaling pathway	5	0.010473786
B2 t16 vs. B2 t24	Down	GO:0002540	leukotriene production involved in inflammatory response	2	0.002222592
B2 t16 vs. B2 t24	Down	GO:0019370	leukotriene biosynthetic process	2	0.007759634
B2 t16 vs. B2 t24	Up	GO:0060612	adipose tissue development	5	2.03E-04
B2 t16 vs. B2 t24	Up	GO:0070373	negative regulation of ERK1 and ERK2 cascade	5	0.003254965
B2 t16 vs. B2 t24	Up	GO:0007264	small GTPase mediated signal transduction	10	0.00727077
B2 t16 vs. B2 t24	Up	GO:0031589	cell-substrate adhesion	3	0.016133553
B2 t16 vs. B2 t24	Up	GO:0033138	positive regulation of peptidyl-serine phosphorylation	5	0.016364796
B2 t16 vs. B2 t24	Up	GO:0016477	cell migration	7	0.01765893
B2 t16 vs. B2 t24	Up	GO:0043277	apoptotic cell clearance	3	0.019873166
B2 t16 vs. B2 t24	Up	GO:0032720	negative regulation of tumor necrosis factor production	3	0.043108056

Enrichment analysis was performed with gene expression data associated with the B2 haplotype. Enrichment was calculated for genes exhibiting statistically significant differences in expression across the successive time points (S1 Table). Three distinct annotation databases were used for enrichment analysis: Gene Ontology—Biological Process, KEGG Pathways, and Reactome Pathways. Complete enrichment annotation available in S2, S3 and S4 Tables.

Between the -6 day and -3 day time points, a large number of genes exhibit reduced expression. These genes are enriched for biological processes such as transcription, mRNA splicing, tRNA processing and negative regulation NFκB mediated gene expression. Similarly, KEGG pathways enriched include RNA polymerase, RNA degradation, Metabolic pathways and Basal transcription factors. The reactome enriched pathways mirror these results with annotations of RNA polymerase II initiation, mRNA capping, and some immune functions including antigen presentation and cytosolic sensors of pathogen-associated DNA.

The transition from day -3 to t0 (just prior to IFNγ stimulation) correlates with down regulation of genes associated with chromosome segregation, mitotic nuclear division and DNA repair, cell cycle pathways, acetylation and some immunological functions of interleukin 3 and 5 signaling and interleukin receptor signaling. Genes exhibiting increases in expression during this interval were enriched in processes and pathways related to interleukin 4, biosynthesis of amino acids, and metabolic pathways. Within an hour of IFNγ stimulation genes exhibiting increased expression were associated with inflammatory and defense responses, toll-like receptor signaling pathways, cell chemotaxis, MyD88 signaling, Jak-STAT signaling, cell adhesion, FoxO signaling, and raf activation.

Genes exhibiting an increase in expression within two hours of IFNγ stimulation are enriched for biological processes of intracellular protein transport, ER-to-Golgi vesical mediated transport, endosome to lysosomal transport, chromatin remodeling, histone H3 acetylation, regulation of vesicle fusion and protein import into the nucleus. Among the KEGG pathways that exhibit enrichment for these genes are RNA transport, protein export, lysosome, mRNA surveillance pathway, endocytosis and mTOR signaling. Similarly, the enriched reactome pathways mirror these processes and include cytosolic sensors of pathogen-associated DNA, RNA polymerase II initiation and promoter clearance, mTORC1-mediated signaling, MAP2K and MAPK activation, downstream TCR signaling, FCε Receptor 1 mediated MAPK activation, and Clec7A signaling.

Genes exhibiting decreased expression between 2 hours and 4 hours following IFNγ stimulation correspond to reduced expression of cell-cycle pathways and mediators, as well as genes implicated in G1/S transcription, DNA replication, and separation of sister chromatids. Biological processes identified within these genes include mitotic spindle assembly, chromosome segregation, and mitotic spindle checkpoint assembly. Conversely, genes exhibiting increased expression during this same time are enriched for biological processes of inflammatory response, defense response to virus, positive regulation of interleukin 12 production, negative replication of viral genome replication, bacterial defense processes and positive regulation of IL1α secretion. KEGG pathways associated with these genes include cytokine-cytokine receptor interaction and genes implicated in influenza A signaling. Reactome processes identified included stem cell population maintenance and TOR signaling.

Over the remaining time points, from 4 hours to 8 hours, from 8 hours to 16 hours and from 16 hours to 24 hours the B2 cells exhibit a systemic down regulation of the genes that were initially activated during the IFNγ stimulation. Overall, the gene enrichment analysis of the RNA sequence data provides a cellular-level picture of the specific biological processes that occur over time following activation of monocyte-derived macrophages.

Discussion

Previous work in our laboratories investigated the association between chicken haplotype and disease resistance, specifically the enhanced resistance of B2 haplotypes to avian coronavirus IBV [10] and the influence of innate immunity leading to decreased clinical signs of illness. We showed that macrophages play an important role in this enhanced immunity, demonstrating much better activation in response to stimulation [13]. To analyze the gene expression involved in this process leading to increased macrophage nitric oxide release in B2 haplotypes, we stimulated macrophages from B2 and B19 chicks for RNA sequencing. In addition, we had observed different cell morphology when isolated monocytes from B2 and B19 haplotypes were differentiating into macrophages and therefore time points before stimulation, during differentiation of the macrophages, were included in this study.

The rationale for investigating the gene expression differences between the B2 and B19 haplotype birds was to address underlying questions that were raised at the end of our previous studies: 1. Why do the IFNγ stimulated B19 derived macrophages exhibit decreased nitric oxide production compared to the IFNγ stimulated B2 derived macrophages? 2. How do the two lineages of macrophages differ at the gene expression level? 3. What specific patterns of gene expression correlate with divergent macrophage differentiation, activation and function? 4. What is the underlying cause of the divergent gene expression patterns observed between B2 and B19 macrophages?

The data collection, analysis and interpretation of results described herein provide plausible answers to these questions based on bioinformatics and functional genomics approaches. Although these answers are more realistically new hypotheses for further investigation, they do represent significant advances in the understanding of B2 and B19 monocyte differentiation into macrophages and the resulting divergent patterns of B2 and B19 macrophage activation and function.

Ultimately, the findings and interpretations we report must be functionally and experimentally validated. Even so, the use of computational methods to answer these questions represents a valuable first step in deciphering the cellular phenotypes underlying MHC haplotype variation in macrophage cells.

Our results demonstrate that there are large numbers of genes differentially expressed in the two haplotypes, both during differentiation of peripheral monocytes into mature macrophages, as well as after stimulation of differentiated macrophages with interferon.

The answer to the question of why the IFNγ stimulated B2 haplotype cells produce more nitric oxide than the IFNγ stimulated B19 cells lies in the timing of macrophage differentiation and the phenotypic variation that is set up early prior to IFNγ stimulation, such as divergent expression of genes involved in differentiation and immune competence. At day t-6, after plating of monocytes without interferon stimulation, several genes relating to inflammation, interferon responses and differentiation are upregulated in the B2 haplotype. This is to be expected as adherence of the monocytes is actually an activation signal, but it is notable that this signal is not resulting in the same gene expression pattern in the B19 haplotype. This pattern can be observed for genes such as IL1β, PTSG2, IL6, which are mainly associated with the inflammatory M1 phenotype. On the other hand, Adenosine receptor A2B is also showing increased gene expression at this time in B2 cells, and this receptor plays an important role in differentiation, as well as in the inflammatory response. Expressions of these genes are consequently initially increased at the time of adherence, and then again after stimulation with interferon, in the B2 haplotype. Some, but not all of these genes are expressed after stimulation with interferon in the B19 haplotype, but not to the same extent as the B2 macrophages, which appears to relate to the initial lack of expression at t-6 days, the beginning of differentiation.

Another interesting observation was the differential expression of genes at day t-6 versus day t-3 in the two haplotypes, as a large number of genes is highly expressed in both haplotypes at day t-6, but then is completely shut down in B2 haplotypes while showing delayed expression until day t-3 in the B19 birds. This seems to indicate a lack of appropriate regulation in the B19 birds, consequently leading to less coherent initiation of gene expression when stimulated. Some of the genes showing this pattern are macrophage differentiation associated GATA2 and FADD, as well as macrophage podosome markers VCL and GSN. Taken together, these results appear to suggest that the regulation of B2 differentiation from monocyte to macrophage is very tightly regulated with many genes increasing in expression, quickly followed by shutting down this increased gene expression. In contrast, the regulation of B19 gene expression is not well regulated, appearing to “linger” with either delayed or extended gene expression. Consequently, we observed differences in expression of genes after stimulation with interferon. Specifically, genes that were strongly expressed at day t-6 and not expressed (or only weakly expressed) at day t-3 in the B2 birds, were robustly increased at 2 and 4 hours of stimulation. In contrast, the same genes showed weak and delayed expression in B19 birds after stimulation, emphasizing the importance of the regulation of gene expression during differentiation. This relates very well to the differences we previously reported in morphology of B2 and B19 macrophages during differentiation and after stimulation.

Our results provide insight into the complexity associated with macrophage differentiation, activation and function. Thousands of genes are up-regulated and then down-regulated in a 24-hour period following IFNγ stimulation. The coordinated activity of multiple regulatory and gene expression control mechanisms is required to effectively achieve the dramatic changes in internal cellular programing that occur. Although the B2 and B19 birds’ haplotypes differ within the MHC locus, the functional consequences of this genetic difference extend well beyond the genes encoded within the B-locus, including macrophage differentiation, M1 and M2 macrophage markers, lysosomal factors involved in phagocytosis, podosome development, invadosome capabilities, chemotaxis potential and matrix degradation ability. Additionally, during the process of differentiation and activation, thousands of genes associated with basic cellular biology undergo rapid changes in expression in coordination with the expression of factors associated with cell renewal and proliferation such as cell cycle regulators, mitotic spindle components, factors involved in chromatin remodeling, molecules required for chromosome segregation and nuclear division.

Taken together, our data and interpretations provide a framework of possible mechanisms of B2 and B19 macrophage biology in differentiation and activation. As such, our findings offer a number of hypotheses about macrophage cell biology that can be used for subsequent studies aimed at validating our findings. Although we performed RT-PCR on a set of differentially expressed genes between the B2 and B19 macrophages, it is not feasible, or possible, to systematically verify, via RT-PCR, each and every transcript, observed at each time point, in the experiment. Even so, our PCR validation provides independent evidence that the pattern of gene expression we observed in the RNA sequence data, was consistent and reproducible which is also in line with our previous research detailing differences in macrophage activation and function [13]

Conclusions

We have tried to elucidate possible mechanisms involved in enhanced disease resistance and macrophage functions displayed by B2 haplotype chickens compared to B19. This study highlights the complex gene expression patterns involved in macrophage differentiation and activation.

One of the main conclusions from the large number of differences seen in the gene expression of the two haplotypes is the fact that there are not just a few genes or genetic markers that can be readily identified as being the ultimate cause of enhanced macrophage function in B2 chickens. Rather, it appears that events during differentiation of monocytes into macrophages have a significant impact on the subsequent ability for stimulation of immune genes after IFNγ treatment. The differences in gene expression correlate with the previously observed differences in morphology of the two haplotypes, with B2 macrophages having a more typical macrophage appearance [13].

Considering the global temporal dysregulation of many genes in B19 haplotypes compared to the more resistant B2 chicks, it seems likely that a variety of genomic regulatory mechanisms (such as transcription factors, miRNAs, snoRNAs, ubiquitin mediated proteasomal degradation, and epigenetic regulation) might play a major role in this process which will be further detailed in a future publication. It will be of great interest to further elucidate these mechanisms and their connection to enhanced immunity. Ultimately, our detailed model of macrophage differentiation, activation and function following IFNγ stimulation provides a high resolution molecular map of cellular biology which can be leveraged by other investigators to further explore the role of these genes in immunology.

Significant differences in gene expression with P-Values.

Pairwise gene expression differences between samples (B2 and B19 haplotype) and timepoints (-6 days, -3 days, 0 days, 1 hr, 2 hrs, 4 hrs, 8 hrs, 16 hrs, 24 hrs) are provided with p-values in excel format. The analysis described within this manuscript focused on differences between [1] matched time points between B2 and B19 haplotype chickens (such as B2 1 hr versus B19 1 hr) as well as [2] progressive timepoints within the same haplotype group (such as B2 1hr versus B2 2hr and B19 4 hr versus B19 8 hr). Subsequently the data contained in this supplemental file also focuses predominantly on those comparisons. This file contains a total of 163,043 rows including the header line containing field names (—identifier for each gene either as gene symbol or ensemble geneId; —chromosome number and the start-end base pair location of the gene; and —the “paired” samples compared for significant gene expression, note ‘AB’ corresponds to B2 haplotype and ‘EC’ corresponds to B19 haplotype; —indication that the analysis method performed by cuffdiff program within the cufflinks package was ‘OK’; and —the fpkm values for sample 1 and sample 2 respectively; —the log of the ratio of fpkm1 and fpkm2; testStat—the test statistic generated during the statistical analysis; —the p-value corresponding to the difference in expression between sample1 and sample2; —a multiple testing corrected p-value; —‘yes’ indicates that the pairwise difference in expression is in fact statistically significant).

(XLSX)

Click here for additional data file.

Significant gene ontology biological process enrichment analysis results.

Enriched gene ontology biological process terms identified within up and down regulated genes within the B2 haplotype chicken samples across progressive time points (-6 days, -3 days, 0 days, 1 hr, 2 hrs, 4 hrs, 8 hrs, 16 hrs, 24 hrs) are provided in excel file format. Since the B2 chickens exhibited the most robust macrophage phenotype these samples were used for the analysis as a means of characterizing the biological processes that were associated with the altered gene expression across the experimental time points. Note ‘AB’ indicates B2 haplotype. The file contains a total of 362 rows including the header line containing field names (—indicates the specific pair of time points for which gene expression changes were identified; —indicates the specific set of differentially expressed genes, ‘Down’ or ‘Up’; —indicates the specific subset of terms that were used for the analysis; —provides the gene ontology identifier for the identified gene ontology term/annotation; —the specific enriched gene ontology biological process term; —the number of genes within the differentially expressed genes that are mapped to the particular enriched gene ontology term; —the corresponding percent associated with the specific number of genes; —the p-value associated with the gene ontology term enrichment).

(XLSX)

Click here for additional data file.

Significant KEGG pathway enrichment analysis results.

Since the B2 chickens exhibited the most robust macrophage phenotype these samples were used for the analysis as a means of characterizing the KEGG pathways that were associated with the altered gene expression across the experimental time points. Note ‘AB’ indicates B2 haplotype. The file contains a total of 110 rows including the header line containing field names (—indicates the specific pair of time points for which gene expression changes were identified; —indicates the specific set of differentially expressed genes, ‘Down’ or ‘Up’; —indicates the specific subset of terms that were used for the analysis; —provides the KEGG pathway identifier for the identified pathway term/annotation; —the specific enriched KEGG pathway term; —the number of genes within the differentially expressed genes that are mapped to the particular enriched KEGG pathway term; —the corresponding percent associated with the specific number of genes; —the p-value associated with the KEGG pathway term enrichment).

(XLSX)

Click here for additional data file.

Significant reactome pathway enrichment analysis results.

Since the B2 chickens exhibited the most robust macrophage phenotype these samples were used for the analysis as a means of characterizing the Reactome pathways that were associated with the altered gene expression across the experimental time points. Note ‘AB’ indicates B2 haplotype. The file contains a total of 110 rows including the header line containing field names (—indicates the specific pair of time points for which gene expression changes were identified; —indicates the specific set of differentially expressed genes, ‘Down’ or ‘Up’; —indicates the specific subset of terms that were used for the analysis; —provides the Reactome pathway identifier for the identified pathway term/annotation; —the specific enriched Reactome pathway term; —the number of genes within the differentially expressed genes that are mapped to the particular enriched Reactome pathway term; —the corresponding percent associated with the specific number of genes; —the p-value associated with the Reactome pathway term enrichment).

(XLSX)

Click here for additional data file.