Marek’s disease (MD) is a neoplastic disease in chickens caused by the MD virus (MDV). Successful vaccine development against MD has resulted in increased virulence of MDV and the understanding of genetic resistance to the disease is, therefore, crucial to long-term control strategies. Also, epigenetic factors are believed to be one of the major determinants of disease response.
Results
Here, we carried out comprehensive analyses of the epigenetic landscape induced by MDV, utilizing genome-wide histone H3 lysine 4 and lysine 27 trimethylation maps from chicken lines with varying resistance to MD. Differential chromatin marks were observed on genes previously implicated in the disease such as MX1 and CTLA-4 and also on genes reported in other cancers including IGF2BP1 and GAL. We detected bivalent domains on immune-related transcriptional regulators BCL6, CITED2 and EGR1, which underwent dynamic changes in both lines as a result of MDV infection. In addition, putative roles for GAL in the mechanism of MD progression were revealed.
Conclusion
Our results confirm the presence of widespread epigenetic differences induced by MD in chicken lines with different levels of genetic resistance. A majority of observed epigenetic changes were indicative of increased levels of viral infection in the susceptible line symptomatic of lowered immunocompetence in these birds caused by early cytolytic infection. The GAL system that has known anti-proliferative effects in other cancers is also revealed to be potentially involved in MD progression. Our study provides further insight into the mechanisms of MD progression while revealing a complex landscape of epigenetic regulatory mechanisms that varies depending on host factors.
Rapid advances in epigenetics have led to the discovery of complex mechanisms of gene regulation involving phenomena such as DNA methylation and chromatin modifications. Methylation of particular histone residues has been found to correlate with specific and often opposing cellular functions, e.g. trimethylation of histone H3 lysine 4 (H3K4me3) is associated with transcriptional start sites (TSSs) of active genes while trimethylation of histone H3 lysine 27 (H3K27me3) is found to mark broad genomic regions for repression. Recent studies have also suggested that characteristic combinations of histone modifications or ‘chromatin states’ define functional elements of the genome and determine their contribution to transcriptional regulation[1–3]. Moreover, the epigenetic state of host genes have been observed to be affected by viral infection leading to tumors in humans[4–6]. Thus, epigenetics constitute a dynamic regulatory framework linking genotypes with environmental factors that could play a major role in differential disease responses among individuals having high genetic similarity.
Marek’s disease (MD) is a highly contagious disease caused by an oncogenic α-herpesvirus MD virus (MDV) and characterized by T-cell lymphomas in chickens[7]. Major losses to the poultry industry as a result of MD have largely been averted due to the success of various vaccination strategies which, remarkably, is also the first instance of the successful control of a natural cancer-causing agent using vaccines[8–10]. However, the virulence of the virus may have progressively increased as a consequence of vaccine development[11–13]. Several reported instances of vaccine breaks or reduced efficacy of vaccination, therefore, underlines the importance of investigating resistance to the disease as a long-term strategy to control MDV.
Natural resistance to MDV can be divided into two categories: major histocompatibility complex (MHC)-associated resistance, wherein different MHC haplotypes at the B blood group locus confer varying levels of resistance and non-MHC associated resistance in which birds having the same MHC haplotype exhibit vastly different responses to MDV infection. Inbred lines 63 and 72 developed at the Avian Disease and Oncology Laboratory (ADOL, East Lansing, MI) that we used in this study, fall into the latter category. These lines share a high degree of genetic similarity but have divergent responses to MDV infection completely independent of the MHC. Several studies have attempted to pinpoint factors responsible for conferring resistance[14–16], but confounding factors, such as, tissue types, virus strains and ages of birds have made it difficult to find a consensus. Multiple risk elements are possibly at play in this complex disease, and increased resistance or susceptibility is likely to be produced by a combination of such factors. In this study, we take a closer look at epigenetic factors behind different responses to MD with a view to a deeper understanding of the broader genomic impact of MDV infection.
We utilized the above population of inbred chickens – line 63 is highly resistant to MD, while line 72 is highly susceptible – and compared the epigenetic effects of MD. Genome-wide maps of H3K4me3 and H3K27me3 in thymus tissues of birds from these chicken lines at the latent stage of MDV infection were generated. We carried out systematic analyses to find differential chromatin marks induced by MDV infection. We also investigated co-localization patterns of the two chromatin modifications to detect putative bivalent domains and the effect of MDV on such domains. The results of our study confirm that Marek’s disease has far-reaching effects on the epigenetic landscape of chicken lines with diverse responses to the virus and, thus, furthers our understanding of this complex disease.
Results
Genome-wide distribution of H3K4me3 and H3K27me3
We performed ChIP-Seq experiments on infected and uninfected birds from lines 63 and 72 to investigate the epigenetic effects of MDV infection. More than 76 million reads from eight samples were mapped to the chicken genome yielding 14418 and 24950 significantly enriched regions (SERs) for H3K4me3 and H3K27me3, respectively (Table1). We further classified these regions as follows: Ubiquitous SERs were found in all samples and were likely due to similarities in the genetic background of the chickens. Line-specific SERs were present in only one line either before or after MDV infection, while condition-specific SERs appeared in both lines but only in individuals with the same infection status.
Table 1
Significantly enriched regions (SERs) and associated genes in each sample
H3K4me3
H3K27me3
Samples
SERs (%)
Genes
SERs (%)
Genes
Line-Specific
63I
647 (4.5)
78
3477 (13.9)
615
63N
594 (4.1)
71
2514 (10.1)
896
63I,63N
924 (6.4)
190
577 (2.3)
150
72I
105 (0.7)
16
1658 (6.6)
451
72N
126 (0.9)
11
2506 (10)
346
72I,72N
73 (0.5)
17
330 (1.3)
89
Condition-specific
63I,72I
97 (0.7)
35
2061 (8.3)
579
63N,72N
47 (0.3)
9
66 (0.3)
22
Ubiquitous
63I,63N,72I,72N
10691 (74.2)
9475
5831 (23.4)
2942
Total
14418
10206
24950
7904
63I: line 63 infected, 63N: line 63 control, 72I: line 72 infected, 72N: line 72 control.
Ubiquitous SERs formed the largest percentage of all enriched regions, accounting for 74.2% and 23.3% in H3K4me3 and H3K27me3 samples, respectively. In the case of H3K4me3, there were large differences in the number of specific SERs - more than 2000 line-specific SERs were found in line 63, compared to about 300 in line 72. Similarly, we found 50% more line-specific SERs of H3K27me3 in line 63 (6568) compared to line 72 (4494). However, upon closer examination, most of the line-specific and condition-specific SERs were revealed to have low read counts (Additional file1: Figure S1) corresponding to regions of low enrichment.
Genes were divided into five regions – promoter, 5’ untranslated region (UTR), exons, introns and 3’ UTR – and the distribution of SERs across these elements was probed (Figure1A). We found a large number of intergenic regions marked by H3K27me3, consistent with high levels of this mark associated with areas of silent heterochromatin. In the case of H3K4me3, a larger proportion of SERs were found around the promoter, exons and the 5’ UTR while similar proportions of H3K4me3 and H3K27me3 SERs were present in introns and 3’ UTRs. A comparison of the genomic distributions of SERs in the different samples (Additional file1: Figures S2A, B) showed a similar number of H3K4me3 SERs across the promoter, exons and the 5’ and 3’ UTRs of genes. Line 63 contained a higher number of intronic and intergenic SERs as compared to line 72 although this did not appear to change as a result of MDV infection. On the other hand, a greater number of H3K27me3 SERs were found in the infected samples although these levels were similar in the two different lines.
Figure 1
Genomic distribution of SERs and relationship between histone marks and gene expression. (A) Distribution of SERs over different genomic elements shows increased levels of H3K4me3 around the promoter region and exons while there are increased levels of H3K27me3 on intergenic regions. (B-E) Relationship between gene expression and histone marks in infected line 63 birds. Plots of histone modifications around the gene body (B, C) in genes having high (blue), medium (red), low (green) and no activity (brown): H3K4me3 shows positive correlation with gene expression levels while H3K27me3 exhibits a negative relationship. A comparison of epigenetic marks and transcriptional levels (D, E) confirms the same. Similar trends were observed in other experimental groups (Additional file1: Figures S3-5).
To analyze the relationship between histone modifications and gene expression, histone modification profiles surrounding the TSS and gene body were plotted for genes ranked by their expression level (Additional file1: Figures1B-E and S3-5). As expected, a strong positive correlation was observed between gene expression and H3K4me3 marks with a distinct peak around the TSS and little to no enrichment within the gene body. On the other hand, H3K27me3 showed negative correlation with gene expression with a peak near the TSS followed by a broad plateau across the gene body. However, the latter relationship was non-linear – genes with lower expression had similar levels of H3K27me3 marks and levels were markedly distinct only at higher expression levels (Figure1C, E).
Differential H3K4me3 marks on genes related to MD
We conducted a comprehensive analysis of genome-wide chromatin marks to find significant differences in MDV-induced responses in line 63 and 72. We used two sets of comparisons: First, to assess the influence of MDV infection within each line, we compared the infected and the non-infected samples from the same line. Secondly, the non-infected samples from the two lines were compared to detect line-specific effects. As a result of this analysis we found 179 differential H3K4me3 SERs and 1116 differential H3K27me3 SERs that mapped to 59 and 66 genes, respectively (Table2). A majority of differential SERs were found in the comparison between non-infected samples of the two lines (Additional files2 and3) with several observed on genes that have been associated with MDV infection.
Table 2
Differential SERs identified in thymus
H3K4me3
H3K27me3
Comparison
Differential SERs*
Genes
Differential SERs*
Genes
63I vs 63N
9
4
42
1
72I vs 72N
30
13
5
0
63N vs 72N
148
46
1094
65
Total
179
59†
1116
66
*FDR < 0.4. † Some genes are shared between different comparisons. 63I: line 63 infected, 63N: line 63 control, 72I: line 72 infected, 72N: line 72 control.
MX1 is a zinc-finger transcription factor that has antiviral properties against a large number of viruses. MX1 was upregulated after MDV infection[17] although its contribution to MD progression is unknown. MDV infection induced a highly significant increase in H3K4me3 enrichment in the promoter region of MX1 in both lines (line 63: p = 1.28x10-7, line 72: p = 4.26x10-9; Figure2A). We observed a concurrent increase in transcript levels after MDV infection in line 72 (p = 0.0264; Figure2B); MX1 expression in line 63 showed a similar trend (fold change = 38.22, p=0.085) although mRNA levels were much lower.
Figure 2
Genes related to MD show differential H3K4me3 marks.MX1 (A, B) and CTLA-4 (C, D) show increased H3K4me3 marks and higher expression in infected individuals from both lines; MMP2 (E, F) exhibits higher levels of H3K4me3 in susceptible line 72. n = 4; * = significant at p < 0.05; ** = significant at p < 0.01; *** = significant at p < 0.001.
CTLA-4 is a cell surface glycoprotein expressed on CD4+ and CD8+ T lymphocytes that is a powerful negative regulator of T-cell activation[18]. The CTLA4 protein is expressed on T lymphocytes soon after activation and regulates T-cell proliferation, production of IL-2 and also supports the function of Treg cells that suppress immune response[19]. Previous studies have reported an increase in CTLA-4 expression after MDV infection[20]. We found an increase in H3K4me3 enrichment in line 72 as a result of MDV infection (p = 0.0003) and there was a similar trend in line 63 (Figure2C). Consistent with the above, there was a significant increase in transcript levels after MDV infection in line 72 (p = 0.04) and a small increase in line 63 (Figure2D).
MMP2 plays a key role in the degradation of the extra-cellular matrix, and an increase in expression has been associated with increasing tumor cell migration and tumor angiogenesis[21, 22]. MMP2 was upregulated during the neoplastic stage of MD infection in susceptible birds[23] but downregulated in response to MDV infection during early lytic infection in susceptible and resistant chickens[17]. We observed a slight increase in H3K4me3 enrichment after MDV infection in both lines, while line 72 exhibited significantly lower levels than line 63 (p = 0.0016; Figure2E). This was coupled with increased MMP2 expression in line 63 after infection (p = 0.0068) while there was no such change in line 72 (Figure2F).
Genes related to cancers show epigenetic changes in response to MD
We observed differential histone marks on several genes that have been associated with other cancers but not in the context of MDV infection. Insulin-like growth factor 2 binding protein 1 (IGF2BP1) is an RNA-binding factor that regulates the translation of mRNAs produced by certain genes like IGF2 and ACTB. Increased expression of IGF2BP1 has been implicated in the development and progression of cancers of various organs, e.g. lung, brain, breast and colon[24–27]. There was no change in the H3K4me3 enrichment levels induced by MDV infection although a significantly higher level of enrichment was present in line 72 (p = 4.21x10-13; Figure3A). Transcript levels in line 72 were much higher than in line 63, but reduced in response to MDV infection (p = 0.044) (Figure3B).
Figure 3
MD induces epigenetic changes in genes related to various cancers.IGF2BP1 (A, B) shows differential H3K4me3 marks and increased expression in susceptible birds while EAF2 (C, D) and GAL (E, F) have differential H3K27me3 levels on the gene body. n = 4; * = significant at p < 0.05; ** = significant at p < 0.01; *** = significant at p < 0.001.
ELL associated factor 2 (EAF2) is a testosterone regulated apoptosis inducer and tumor suppressor. Inactivation of EAF2 has been shown to lead to tumors in multiple organs[28]. There was a significant increase in H3K27me3 levels after MDV infection in line 63 (p = 0.0414) while among uninfected chickens these levels were markedly higher in line 72 (p = 0.0138; Figure3C). However, EAF2 expression was drastically reduced after MDV infection in line 72 (p=0.0016) but showed only a small decrease in line 63 (Figure3D).
Galanin (GAL) is a neuropeptide that modulates various physiological functions, such as, inhibition of insulin secretion and stimulation of growth hormone secretion. Three galanin receptors are known (GALR1, 2 and 3): the expression of GALR1 has anti-proliferative effects while GALR2 can be both anti- or pro-proliferative in function. Therefore, the GAL system is considered to be a promising candidate for detection and treatment of various cancers[29, 30]. We observed an increase in H3K27me3 levels on GAL after infection in both lines (Figure3E). Also, expression levels were significantly lowered after MDV infection in line 72 (p = 0.00087) while there was a similar trend in line 63 (p = 0.051; Figure3F). Interestingly, GALR1 had both H3K4me3 and H3K27me3 enrichments (Figure4) although GALR2 showed no significant histone marks (Additional file1: Figure S6).
Figure 4
Significant H3K4me3 and H3K27me3 enrichment aroundGALR1. The anti-proliferative GAL receptor GALR1 exhibited both active and repressive marks. There is no change in H3K4me3 levels but a definite increase in H3K27me3 levels after infection in line 72.
Regions of chromatin having both active and repressive marks are said to be bivalent and have been shown to play important roles in development and genetic imprinting[31, 32]. For example, bivalent domains have been shown to mark promoters of genes that are subsequently silenced in tumors by DNA hypermethylation indicating their importance in cancer[33]. A mono-allelic bivalent chromatin domain that controls tissue-specific genomic imprinting at a specific locus was recently found in mice[32]. To investigate the presence of such bivalent chromatin states and the possible effect of MDV infection, we defined bivalent genes as those having H3K4me3 reads (TSS ± 500 bp) greater than 30 reads per million mapped reads (RPM) and H3K27me3 reads (gene body) greater than 2 RPM, respectively (~85th percentile). This filtering process yielded a list of 99 putative bivalent genes (Additional file4).
Functional annotation clustering of the above genes using DAVID[34, 35] revealed significant enrichment of several immune-related functions. These included transcription factor EGR1 which is reported to have tumor suppressor properties, genes involved in lymphocyte activation and differentiation such as BCL6, CD4 and SMAD3 and genes TLR3 and TIRAP that are part of the toll-like receptor signaling pathway. Bivalent domains were also present on a variety of transcription factors with immune-related functions such as CITED2, a transactivator that regulates NF-κB, MYC a transcription factor associated with hematopoetic tumors and RHOB a Ras family homolog that mediates apoptosis in tumor cells after DNA damage. Moreover, all genes involved in the top five functional annotation clusters showed higher chromatin levels in line 72 primarily after MDV infection (Additional file5).
Bivalent domains are altered in response to MD
We further investigated the effect of MD on bivalent chromatin domains observed on BCL6, CITED2, EGR1, CD4 and TLR3 (Additional file1: Figure5 and S7). Interestingly, three of these genes, CITED2, BCL6 and EGR1, showed identical epigenetic and transcriptional signatures.
Figure 5
Bivalent domains on transcriptional regulators are altered by MD. H3K4me3 (orange) and H3K27me3 (green) profiles and associated transcript levels of BCL6 (A, B), EGR1 (C, D) and CITED2 (E, F). In all three cases we observe a slight increase in H3K27me3 induced by MDV infection in line 72 and a concurrent significant decrease in transcript levels while increase in active and repressive marks appear to cancel each other out in line 63.
CITED2 is a member of the p300/CBP co-activator family that has intrinsic histone acetyltransferase activity and plays a major role in regulating and coordinating multiple complex cellular signals to determine the expression level of a gene[36]. B-cell CLL/lymphoma 6 (BCL6) is a zinc finger protein that functions as a transcriptional repressor which was downregulated at 15 dpi in spleen tissues from F1 progeny (15I5 X 71) of MD-susceptible chickens[20]. EGR1 belongs to a group of early response genes induced by a variety of signaling molecules such as growth factors, hormones and neurotransmitters that is believed to play a major role in cell proliferation and apoptosis[37]. Overexpression of EGR1 promotes tumor growth in kidney cells[38] but suppresses growth and transformation in other cell types, e.g. fibroblasts and glioblastoma cells[39].
In each of the above genes, both active and repressive chromatin marks were increased in response to infection in line 63 chickens. However, in line 72, there was a definite increase in H3K27me3 marks but no change in H3K4me3 (Figures5A, C, E). Transcript levels were in agreement with this observation: infected line 72 chickens showed a significant downregulation (CITED2: p=0.0004; BCL6: p=0.0048; EGR1: p=0.0005; Figures5B, D, F), while there were no such changes in line 63.
On the other hand, TLR3 and CD4 showed a slight increase in H3K4me3 marks after MDV infection while there were no appreciable changes in H3K27me3 levels. In keeping with the epigenetic changes, there was a small increase in gene expression in infected birds from both lines (Additional file1: Figure S7).
Discussion
Immune parameters that are known to play a major role in genetic resistance to MDV are correlated with innate immune responses, such as NK cell activity, production of nitric oxide and cytokines, such as, interferons. Recent studies have identified host cytokines such as IL-18 and IFN-γ that contribute to the initiation and continuation of latency[40]. However, cytokine levels can undergo rapid flux in response to infection, and consistent with this, we did not observe any epigenetic changes associated with these genes in the MHC-congenic lines used in our study (Additional file1: Figure S8). This suggests the existence of other extrinsic factors responsible for transcriptional variations between resistant and susceptible chickens at this stage of the disease.
A global comparison of histone modification levels in the two inbred chicken lines yielded some interesting results. As expected, a majority of SERs were found in all the experimental groups, indicating a high level of epigenetic similarity between the lines in addition to inherent genetic similarity. In the case of H3K27me3, the percentage of ubiquitous SERs was relatively low (23.4%), although this was likely due to lower precision of the peak caller for broad chromatin marks. Besides, most of the SERs detected in a subset of samples corresponded to regions of low enrichment, which may also be the reason behind the relatively low number of differential SERs detected in our study.
Genes that have been previously associated with MD and various other cancers showed differential marks that are either MD-induced (MX1, CTLA-4, EAF2 and GAL) or line-specific (IGF2BP1 and MMP2). The increase in H3K4me3 marks around the TSS of MX1, a gene with known antiviral properties, appeared to be correlated with upregulated expression in both lines in response to MDV infection. However, lowered overall mRNA levels in the resistant line suggest additional factors could be involved in the regulation of this gene. Similarly, increased mRNA levels of the lymphocyte surface marker CTLA4 is possibly due to the presence of larger numbers of T cells in MDV infected birds as a result of higher levels of infection. EAF2 functions as an apoptosis inducer in addition to being a tumor suppressor, and therefore, its downregulation could contribute to higher tumor incidence rates in line 72. However, it is not clear why a significant increase in H3K27me3 levels did not have any effect on transcript levels in the resistant line.
IGF2BP1 is believed to act by stabilizing the mRNA of the c-myc oncogene and therefore, the higher expression of this gene and a more stable c-myc gene product might play a role in increasing MD susceptibility in line 72 birds via increased cell proliferation and transformation. The matrix metalloprotease MMP2 is upregulated after infection in the resistant line 63, similar to the previously observed increase at the neoplastic stage of MD. However, mRNA levels were similar in the two lines before infection contrary to the difference in H3K4me3 levels suggesting that additional factors are responsible for regulating this gene.
The correlation between observed differential histone marks and transcript levels was moderate at best. Indeed, differential H3K4me3 marks were strongly predictive of gene expression levels but the correlation between H3K27me3 and mRNA levels was relatively poor. There could be various reasons for this – indeed, H3K27me3 levels had a non-linear relationship with gene expression with higher marks showing little difference in the effect on expression. Therefore, in this tissue, the levels of H3K27me3 may not be a very good indicator of gene expression levels. Also, the transcription of these genes might be controlled by other factors with the change in H3K27me3 levels only incidental.
Bivalent domains were detected on transcriptional regulators BCL6, CITED2 and EGR1 and the galanin receptor GALR1. The epigenetic and transcriptional signatures observed on these genes indicated that they were poised at the latent stage of the disease, but with crucial differences in the two lines. Increased repressive marks in the susceptible line correlated with significant downregulation of the genes, while in line 63, the increase in both marks appeared to compensate for each other with no eventual effect on gene transcription. This suggests that some ‘poised’ bivalent genes can become highly repressed even with a relatively small increase in H3K27me3 marks. The change in the chromatin levels could also be correlated with an increase in cell populations having the repressive mark. Taken together, the above findings point towards the existence of finely balanced epigenetic control of transcription, which may be necessary to mount a rapid response by the immune system. However, this machinery could potentially be hijacked by a pathogen and result in an aberrant phenotype. The effect of MDV infection on the bivalent domain on GALR1, in particular, suggests the repression and potential loss of its anti-proliferative effects. Thus, the galanin system possibly plays an important and hitherto unknown role in MD progression and could be a novel target for long-term control of the disease.
One of the major players in MDV-induced malignant transformation is Meq, a virus-encoded oncoprotein that has diverse functions, e.g. transactivation, chromatin remodeling and regulation of transcription. Meq interacts with and sequesters the tumor suppressor protein p53, resulting in the dysregulation of cell-cycle control[7] and inhibition of the transcriptional and apoptotic activities of the protein[41]. Several genes that show epigenetic changes in response to MDV infection have been associated with p53. Downregulation of CITED2 stabilizes the p53 protein leading to its accumulation[42]. The BCL6 gene product is believed to contribute to lymphomagenesis by inactivation of p53[43]. Besides, EAF2 has also been shown to interact with and suppress the function of p53[28]. The downregulation of all of the above genes in susceptible birds after MDV infection points towards the increased production of p53 and a robust anti-tumor response. That we still observe higher tumor incidence rates in this line, suggests the presence of large amounts of inactivating viral Meq protein which, in turn, indicates that increased numbers of MD-infected cells are present in the susceptible line at this stage of the disease. A majority of tumors are believed to result from the viral transformation of CD4+ T cells some of which are infected at the latent stage of MD[44]. The larger number of virus-infected cells produced in the susceptible line is possibly due to lowered immunocompetence as a result of the early stages of infection. Thus, a more detailed investigation of the early cytolytic stage of MD is necessary to shed further light on the causes behind the divergent response to MD observed in these birds.
Whole tissues represent a mixture of various cell populations, and observed epigenetic changes might be due to a change in chromatin marks in a particular cell type or a variation in the relative number of cells carrying active or repressive histone marks. However, such in vivo studies are representative of the host response at a systems level wherein different cell types might interact in a cooperative manner to fight infection. Thus, while the study of pure cell populations is likely to yield greater discriminative power, the investigation of tissue macroenvironments is, perhaps, closer to reality.
This study focused on the thymus tissue as it is a major immune organ and contains a significant population of T lymphocytes in various stages of differentiation. Our earlier study of the MDV-induced transcriptome in these birds indicated the presence of line-specific differences at the latent stage of MD[45]. In addition, birds susceptible to MD suffer thymic atrophy during the early stages of infection[46], indicating the importance of understanding changes in this organ to elucidate the mechanisms involved in disease progression. Ongoing studies in our lab include other tissues, e.g. spleen[47], and a time-course through the various stages of infection, to further investigate the systemic effects of MD and the epigenetic basis of MD resistance.
Conclusions
We studied the effect of latent MDV infection on two chromatin modifications in inbred chicken lines exhibiting different degrees of resistance to MD. Several genes showed changes in histone enrichment and this response was often significantly different between the two chicken lines. A detailed analysis of co-localization patterns of the chromatin marks revealed the presence of bivalent domains on a number of immune-related transcriptional regulators. More importantly, these domains showed marked changes in response to MDV infection and provide further evidence of the far-reaching epigenetic effects of MD. Our results suggest putative roles for the GAL system in MD progression. A majority of the differential chromatin marks are also suggestive of increased levels of viral infection in the susceptible line symptomatic of lowered immunocompetence in these birds at early stages of the disease. In summary, our study provides further insight into the mechanisms of MD progression while revealing a complex landscape of epigenetic regulatory mechanisms. Further work is necessary to fully elucidate the underlying mechanisms of MD, but our results suggest that this is a promising step towards a deeper understanding of this complex disease.
Methods
Animals and viruses
Two specific-pathogen-free inbred lines of White Leghorn either resistant (63) or susceptible (72) to MD were hatched, reared and maintained in Avian Disease and Oncology Laboratory (ADOL, Michigan, USDA). Four chickens from each line were injected intra-abdominally with a partially attenuated very virulent plus strain of MDV (648A passage 40) at 5 days after hatch with a viral dosage of 500 plaque-forming units (PFU). Infected and control chickens from both lines (n = 4) were terminated at 10dpi to collect thymus tissues. All procedures followed the standard animal ethics and use guidelines of ADOL.
Quantification of MDV loads in thymus
The MDV gene ICP4 was used for quantification of viral genomic DNA in thymus as previously described[48]. Quantitative PCR was performed by using 100 ng/μl of genomic DNA on the iCycler iQ PCR system (Bio-Rad, USA) and QuantiTect SYBR Green PCR Kit (Qiagen, USA) (Additional file1: Figure S9). The relative MDV loads were determined by normalizing to a single-copy gene Vim (vimentin)[49]. The primers for Vim are as follows: Forward: 5’-CAGCCACAGAGTAGGGTAGTC-3’; Reverse: 5’-GAATAGGGAAGAACAGGAAAT-3’.
Chromatin Immunoprecipitation and Illumina Sequencing
Chromatin immunoprecipitation (ChIP) was carried out using thymus samples from MDV infected and controls birds[50]. About 30 mg thymus samples from three individuals were cut into small pieces (1 mm3) and digested with MNase to obtain mononucleosomes. PNK and Klenow enzymes (NBE, Ipswich, MA, USA) were used to repair the ChIP DNA ends pulled down by the specific antibody. A 3′ adenine was then added using Taq polymerase and a pair of Solexa adaptors (Illumina, USA) ligated to the repaired ends. Seventeen cycles of PCR was performed on ChIP DNA using the adaptor primers and fragments with a length of about 190 bp (mononucleosome + adaptors) were isolated from agarose gel. Subsequently, cluster generation and ChIP-sequencing (ChIP-Seq) using the purified DNA was performed on the Solexa 1G Genome Analyzer (Illumina, USA) following manufacturer protocols. The antibodies used and the total number of reads obtained for each sample is listed in Additional file6.
Read mapping and summary counts
Sequence reads obtained from the Illumina 1G Genome Analyzer were aligned to the May 2006 version of the chicken genome (galGal3) using Maq version 0.7.1[51]. Default alignment policies of Maq were enforced: a valid alignment could have a maximum of two mismatches and if a read aligned equally well to multiple places in the genome, one was chosen at random. If multiple reads mapped to the same genomic location only one was kept to avoid amplification bias. Summary read counts were calculated using non-overlapping windows of 200 bp for visualization and normalized to per million mapped reads in each sample for the purpose of comparisons.
Identification of enriched regions
Summarized read counts were subjected to peak calling with SICER[52]. The source code was modified to include support for the chicken genome. Fragment length was specified to be 190. A window size of 200 bp and gap size of 400 bp was used for the analysis. The E-value for estimating significant peaks was set to 100. For the purposes of comparing different samples, SERs found in similar genomic regions of different samples were merged to obtain a consolidated list as follows: SERs from one sample were used to initialize the list. For each such region M, we searched for overlapping SERs in the next sample. In the case of an overlap between M and a significant region, S, the merged region was updated to include both M and S. This procedure was iterated over all samples to obtain a consolidated list of merged SERs.
Gene annotation and genomic distribution of SERs
RefSeq and Ensembl gene annotations were downloaded from UCSC genome browser[53]. As there were only 4306 RefSeq genes in the database, we included Ensembl genes in our analysis to improve genome-wide coverage. There were 17858 annotated genes in the Ensembl database, which include validated and predicted genes. Redundancies between the databases were listed and accounted for, yielding a non-redundant list of 18198 genes with 4306 RefSeq genes and 13892 Ensembl genes. We then searched for overlaps between merged SERs and the non-redundant list of annotated genes. For H3K4me3, an SER was annotated with a gene if it overlapped the TSS region of the gene whereas in the case of H3K27me3, a valid overlap constituted an SER overlapping the gene body. To calculate the genomic distribution we counted all SERs having an overlap with one of the following regions: promoter (TSS ± 1 kb), exons, introns, 5’ UTR and 3’ UTR.
Histone modification profiles and differential chromatin marks
Genes were divided into 10 sets based on their absolute expression and representative sets corresponding to high, medium, low and no expression were chosen for visualization. We defined the gene body as the region between the transcription start site (TSS) and the transcription termination site (TTS). Histone modification profiles surrounding the gene body were calculated in 3 distinct regions: 5000 bp upstream of the 5’ end, gene body and 5000 bp downstream of the 3’ end of the gene. For reads falling within the gene body, read counts were obtained in bins 5% of the gene length while 1000 bp windows were used for the 5’ and 3’ flanking regions. The read counts in all cases were normalized to the total number of genes in the categories and total number of reads in the sample. We also compared gene expression to histone modification levels by plotting normalized microarray data (Zhang, H. unpublished data) against reads mapping to (i) TSS ± 500 bp and (ii) the gene body for each gene.
Reads mapping to merged SERs were tested for epigenetic changes induced by MDV infection in lines 63 and 72 using DESeq[54]. We used the method ‘blind’ for variance estimation and p-values were corrected for multiple testing using the Benjamini-Hochberg FDR procedure[55]. Statistical significance was defined using a false discovery rate (FDR) threshold of 0.4.
Validation of ChIP, ChIP-Seq and gene transcription by Q-PCR
Quantitative real-time RT-PCR was used to validate the quality of the ChIP and gene transcript levels on the iCycler iQ PCR system (Bio-Rad, Hercules, CA, USA). The real-time RT-PCR reactions were performed with a QuantiTect SYBR Green PCR Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. An online primer system (http://frodo.wi.mit.edu/primer3/) was used to design the primers and four biological and four technical replicates were performed. The primer sequences are shown in Additional file7.
Data access
Raw and processed ChIP-Seq data discussed in this manuscript were deposited in the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) under Series accession number GSE24017.
Author contributions
AM performed the data analysis and wrote the manuscript. JL performed the ChIP experiments, sequence library preparation, validation of ChIP-Seq results and revised the manuscript. HMZ collected samples and revised the manuscript. JZS designed the experiments and revised the manuscript. All authors read and approved the final manuscript.
Notes
Declarations
Acknowledgements
This project was supported by National Research Initiative Competitive Grant no. USDA-NRI/NIFA 2008-35204-04660 and USDA-NRI/NIFA 2010-65205-20588 from the USDA National Institute of Food and Agriculture.
Authors’ Affiliations
(1)
Department of Animal & Avian Sciences, University of Maryland
(2)
USDA, ARS, Avian Disease and Oncology Laboratory
(3)
Laboratory of Molecular Immunology, National Heart, Lung and Blood Institute, National Institutes of Health
References
Ernst J, Kellis M: Discovery and characterization of chromatin states for systematic annotation of the human genome. Nat Biotechnol. 2010, 28 (8): 817-825. 10.1038/nbt.1662.PubMed CentralView ArticlePubMedGoogle Scholar
Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, Zhang X, Wang L, Issner R, Coyne M, et al: Mapping and analysis of chromatin state dynamics in nine human cell types. Nature. 2011, 473 (7345): 43-49. 10.1038/nature09906.PubMed CentralView ArticlePubMedGoogle Scholar
Weishaupt H, Sigvardsson M, Attema JL: Epigenetic chromatin states uniquely define the developmental plasticity of murine hematopoietic stem cells. Blood. 2010, 115 (2): 247-256. 10.1182/blood-2009-07-235176.View ArticlePubMedGoogle Scholar
Knight JS, Lan K, Subramanian C, Robertson ES: Epstein-Barr virus nuclear antigen 3C recruits histone deacetylase activity and associates with the corepressors mSin3A and NCoR in human B-cell lines. J Virol. 2003, 77 (7): 4261-4272. 10.1128/JVI.77.7.4261-4272.2003.PubMed CentralView ArticlePubMedGoogle Scholar
Tsai CN, Tsai CL, Tse KP, Chang HY, Chang YS: The Epstein-Barr virus oncogene product, latent membrane protein 1, induces the downregulation of E-cadherin gene expression via activation of DNA methyltransferases. Proc Natl Acad Sci U S A. 2002, 99 (15): 10084-10089. 10.1073/pnas.152059399.PubMed CentralView ArticlePubMedGoogle Scholar
Osterrieder N, Kamil JP, Schumacher D, Tischer BK, Trapp S: Marek's disease virus: from miasma to model. Nat Rev Microbiol. 2006, 4 (4): 283-294. 10.1038/nrmicro1382.View ArticlePubMedGoogle Scholar
Epstein MA, Achong BG, Churchill AE, Biggs PM: Structure and development of the herpes-types virus of Marek's disease. J Natl Cancer Inst. 1968, 41 (3): 805-820.PubMedGoogle Scholar
Churchill AE, Payne LN, Chubb RC: Immunization against Marek's disease using a live attenuated virus. Nature. 1969, 221 (5182): 744-747. 10.1038/221744a0.View ArticlePubMedGoogle Scholar
Okazaki W, Purchase HG, Burmester BR: Protection against Marek's disease by vaccination with a herpesvirus of turkeys. Avian Dis. 1970, 14 (2): 413-429. 10.2307/1588488.View ArticlePubMedGoogle Scholar
Davison F, Nair V: Use of Marek's disease vaccines: could they be driving the virus to increasing virulence?. Expert Rev Vaccines. 2005, 4 (1): 77-88. 10.1586/14760584.4.1.77.View ArticlePubMedGoogle Scholar
Gimeno IM: Marek's disease vaccines: a solution for today but a worry for tomorrow?. Vaccine. 2008, 26 (Suppl 3): C31-C41.View ArticlePubMedGoogle Scholar
Witter RL: Increased virulence of Marek's disease virus field isolates. Avian Dis. 1997, 41 (1): 149-163. 10.2307/1592455.View ArticlePubMedGoogle Scholar
Liu HC, Cheng HH, Tirunagaru V, Sofer L, Burnside J: A strategy to identify positional candidate genes conferring Marek's disease resistance by integrating DNA microarrays and genetic mapping. Anim Genet. 2001, 32 (6): 351-359. 10.1046/j.1365-2052.2001.00798.x.View ArticlePubMedGoogle Scholar
Liu HC, Kung HJ, Fulton JE, Morgan RW, Cheng HH: Growth hormone interacts with the Marek's disease virus SORF2 protein and is associated with disease resistance in chicken. Proc Natl Acad Sci U S A. 2001, 98 (16): 9203-9208. 10.1073/pnas.161466898.PubMed CentralView ArticlePubMedGoogle Scholar
Wain HM, Toye AA, Hughes S, Bumstead N: Targeting of marker loci to chicken chromosome 16 by representational difference analysis. Anim Genet. 1998, 29 (6): 446-452. 10.1046/j.1365-2052.1998.296374.x.View ArticlePubMedGoogle Scholar
Smith J, Sadeyen JR, Paton IR, Hocking PM, Salmon N, Fife M, Nair V, Burt DW, Kaiser P: Systems Analysis of Immune Responses in Marek's Disease Virus Infected Chickens Identifies a Gene Involved in Susceptibility and Highlights a Possible Novel Pathogenicity Mechanism. J Virol. 2011, 85 (21): 11146-11158. 10.1128/JVI.05499-11.PubMed CentralView ArticlePubMedGoogle Scholar
McCoy KD, Le Gros G: The role of CTLA-4 in the regulation of T cell immune responses. Immunol Cell Biol. 1999, 77 (1): 1-10. 10.1046/j.1440-1711.1999.00795.x.View ArticlePubMedGoogle Scholar
Curran MA, Montalvo W, Yagita H, Allison JP: PD-1 and CTLA-4 combination blockade expands infiltrating T cells and reduces regulatory T and myeloid cells within B16 melanoma tumors. Proc Natl Acad Sci U S A. 2010, 107 (9): 4275-4280. 10.1073/pnas.0915174107.PubMed CentralView ArticlePubMedGoogle Scholar
Fingleton B: Matrix metalloproteinases: roles in cancer and metastasis. Front Biosci. 2006, 11: 479-491. 10.2741/1811.View ArticlePubMedGoogle Scholar
Rath NC, Parcells MS, Xie H, Santin E: Characterization of a spontaneously transformed chicken mononuclear cell line. Vet Immunol Immunopathol. 2003, 96 (1–2): 93-104.View ArticlePubMedGoogle Scholar
Chen C, Li H, Xie Q, Shang H, Ji J, Bai S, Cao Y, Ma Y, Bi Y: Transcriptional profiling of host gene expression in chicken liver tissues infected with oncogenic Marek's disease virus. J Gen Virol. 2011, 92: 2724-2733. 10.1099/vir.0.034066-0.View ArticlePubMedGoogle Scholar
Ioannidis P, Kottaridi C, Dimitriadis E, Courtis N, Mahaira L, Talieri M, Giannopoulos A, Iliadis K, Papaioannou D, Nasioulas G, et al: Expression of the RNA-binding protein CRD-BP in brain and non-small cell lung tumors. Cancer Lett. 2004, 209 (2): 245-250. 10.1016/j.canlet.2003.12.015.View ArticlePubMedGoogle Scholar
Ioannidis P, Mahaira L, Papadopoulou A, Teixeira MR, Heim S, Andersen JA, Evangelou E, Dafni U, Pandis N, Trangas T: 8q24 Copy number gains and expression of the c-myc mRNA stabilizing protein CRD-BP in primary breast carcinomas. Int J Cancer. 2003, 104 (1): 54-59. 10.1002/ijc.10794.View ArticlePubMedGoogle Scholar
Kato T, Hayama S, Yamabuki T, Ishikawa N, Miyamoto M, Ito T, Tsuchiya E, Kondo S, Nakamura Y, Daigo Y: Increased expression of insulin-like growth factor-II messenger RNA-binding protein 1 is associated with tumor progression in patients with lung cancer. Clin Cancer Res. 2007, 13 (2 Pt 1): 434-442.View ArticlePubMedGoogle Scholar
Ross J, Lemm I, Berberet B: Overexpression of an mRNA-binding protein in human colorectal cancer. Oncogene. 2001, 20 (45): 6544-6550. 10.1038/sj.onc.1204838.View ArticlePubMedGoogle Scholar
Su F, Pascal LE, Xiao W, Wang Z: Tumor suppressor U19/EAF2 regulates thrombospondin-1 expression via p53. Oncogene. 2010, 29 (3): 421-431. 10.1038/onc.2009.326.PubMed CentralView ArticlePubMedGoogle Scholar
Berger A, Santic R, Hauser-Kronberger C, Schilling FH, Kogner P, Ratschek M, Gamper A, Jones N, Sperl W, Kofler B: Galanin and galanin receptors in human cancers. Neuropeptides. 2005, 39 (3): 353-359. 10.1016/j.npep.2004.12.016.View ArticlePubMedGoogle Scholar
Rauch I, Kofler B: The galanin system in cancer. EXS. 2010, 102: 223-241.PubMedGoogle Scholar
Bernstein BE, Mikkelsen TS, Xie X, Kamal M, Huebert DJ, Cuff J, Fry B, Meissner A, Wernig M, Plath K, et al: A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell. 2006, 125 (2): 315-326. 10.1016/j.cell.2006.02.041.View ArticlePubMedGoogle Scholar
Sanz LA, Chamberlain S, Sabourin JC, Henckel A, Magnuson T, Hugnot JP, Feil R, Arnaud P: A mono-allelic bivalent chromatin domain controls tissue-specific imprinting at Grb10. EMBO J. 2008, 27 (19): 2523-2532. 10.1038/emboj.2008.142.PubMed CentralView ArticlePubMedGoogle Scholar
Rodriguez J, Munoz M, Vives L, Frangou CG, Groudine M, Peinado MA: Bivalent domains enforce transcriptional memory of DNA methylated genes in cancer cells. Proc Natl Acad Sci U S A. 2008, 105 (50): 19809-19814. 10.1073/pnas.0810133105.PubMed CentralView ArticlePubMedGoogle Scholar
da Huang W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009, 4 (1): 44-57.View ArticlePubMedGoogle Scholar
da Huang W, Sherman BT, Lempicki RA: Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009, 37 (1): 1-13. 10.1093/nar/gkn923.View ArticlePubMedGoogle Scholar
Chan HM, La Thangue NB: p300/CBP proteins: HATs for transcriptional bridges and scaffolds. J Cell Sci. 2001, 114 (Pt 13): 2363-2373.PubMedGoogle Scholar
Thiel G, Cibelli G: Regulation of life and death by the zinc finger transcription factor Egr-1. J Cell Physiol. 2002, 193 (3): 287-292. 10.1002/jcp.10178.View ArticlePubMedGoogle Scholar
Scharnhorst V, Menke AL, Attema J, Haneveld JK, Riteco N, van Steenbrugge GJ, van der Eb AJ, Jochemsen AG: EGR-1 enhances tumor growth and modulates the effect of the Wilms' tumor 1 gene products on tumorigenicity. Oncogene. 2000, 19 (6): 791-800. 10.1038/sj.onc.1203390.View ArticlePubMedGoogle Scholar
Huang RP, Liu C, Fan Y, Mercola D, Adamson ED: Egr-1 negatively regulates human tumor cell growth via the DNA-binding domain. Cancer Res. 1995, 55 (21): 5054-5062.PubMedGoogle Scholar
Buscaglia C, Calnek BW: Maintenance of Marek's disease herpesvirus latency in vitro by a factor found in conditioned medium. J Gen Virol. 1988, 69 (Pt 11): 2809-2818.View ArticlePubMedGoogle Scholar
Deng X, Li X, Shen Y, Qiu Y, Shi Z, Shao D, Jin Y, Chen H, Ding C, Li L, et al: The Meq oncoprotein of Marek's disease virus interacts with p53 and inhibits its transcriptional and apoptotic activities. Virol J. 2010, 7: 348-10.1186/1743-422X-7-348.PubMed CentralView ArticlePubMedGoogle Scholar
Wu Z-Z, Sun N-K, Chao CCK: Knockdown of CITED2 using short-hairpin RNA sensitizes cancer cells to cisplatin through stabilization of p53 and enhancement of p53-dependent apoptosis. J Cell Physiol. 2011, 226 (9): 2415-2428. 10.1002/jcp.22589.View ArticlePubMedGoogle Scholar
Phan RT, Dalla-Favera R: The BCL6 proto-oncogene suppresses p53 expression in germinal-centre B cells. Nature. 2004, 432 (7017): 635-639. 10.1038/nature03147.View ArticlePubMedGoogle Scholar
Schat KA, Chen CL, Calnek BW, Char D: Transformation of T-lymphocyte subsets by Marek's disease herpesvirus. J Virol. 1991, 65 (3): 1408-1413.PubMed CentralPubMedGoogle Scholar
Yu Y, Luo J, Mitra A, Chang S, Tian F, Zhang H, Yuan P, Zhou H, Song J: Temporal Transcriptome Changes Induced by MDV in Marek's Disease-Resistant and -Susceptible Inbred Chickens. BMC Genomics. 2011, 12 (1): 501-10.1186/1471-2164-12-501.PubMed CentralView ArticlePubMedGoogle Scholar
Morimura T, Ohashi K, Kon Y, Hattori M, Sugimoto C, Onuma M: Apoptosis and CD8-down-regulation in the thymus of chickens infected with Marek's disease virus. Arch Virol. 1996, 141 (11): 2243-2249. 10.1007/BF01718230.View ArticlePubMedGoogle Scholar
Luo J, Mitra A, Tian F, Chang S, Zhang H, Cui K, Yu Y, Zhao K, Song J: Histone Methylation Analysis and Pathway Predictions in Chickens after MDV Infection. PLoS One. 2012, 7 (7): e41849-10.1371/journal.pone.0041849.PubMed CentralView ArticlePubMedGoogle Scholar
Yunis R, Jarosinski KW, Schat KA: Association between rate of viral genome replication and virulence of Marek's disease herpesvirus strains. Virology. 2004, 328 (1): 142-150. 10.1016/j.virol.2004.07.017.View ArticlePubMedGoogle Scholar
Zehner ZE, Paterson BM: Characterization of the chicken vimentin gene: single copy gene producing multiple mRNAs. Proc Natl Acad Sci U S A. 1983, 80 (4): 911-915. 10.1073/pnas.80.4.911.PubMed CentralView ArticlePubMedGoogle Scholar
Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, Wei G, Chepelev I, Zhao K: High-resolution profiling of histone methylations in the human genome. Cell. 2007, 129 (4): 823-837. 10.1016/j.cell.2007.05.009.View ArticlePubMedGoogle Scholar
Li X, Chiang HI, Zhu J, Dowd SE, Zhou H: Characterization of a newly developed chicken 44K Agilent microarray. BMC Genomics. 2008, 9: 60-10.1186/1471-2164-9-60.PubMed CentralView ArticlePubMedGoogle Scholar
Zang C, Schones DE, Zeng C, Cui K, Zhao K, Peng W: A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics. 2009, 25 (15): 1952-1958. 10.1093/bioinformatics/btp340.PubMed CentralView ArticlePubMedGoogle Scholar
Fujita PA, Rhead B, Zweig AS, Hinrichs AS, Karolchik D, Cline MS, Goldman M, Barber GP, Clawson H, Coelho A, et al: he UCSC Genome Browser database: update 2011. Nucleic Acids Res. 2010, 39 (Database issue): D876-D882.PubMed CentralPubMedGoogle Scholar
Benjamini Y, Hochberg Y: Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J Royal Stat Soc Series B (Methodological). 1995, 57 (1): 289-300.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
