Transcriptional Profiling of MEq-Dependent Genes in Marek’s Disease Resistant and Susceptible Inbred Chicken Lines

Marek’s disease (MD) is an economically significant disease in chickens caused by the highly oncogenic Marek’s disease virus (MDV). Understanding the genes and biological pathways that confer MD genetic resistance should lead towards the development of more disease resistant commercial poultry flocks or improved MD vaccines. MDV mEq, a bZIP transcription factor, is largely attributed to viral oncogenicity though only a few host target genes have been described, which has impeded our understanding of MDV-induced tumorigenesis. Given the importance of mEq in MDV-induced pathogenesis, we explored the role of mEq in genetic resistance to MDV. Using global transcriptome analysis and cells from MD resistant or susceptible birds, we compared the response to infection with either wild type MDV or a nononcogenic recombinant lacking mEq. As a result, we identified a number of specific genes and pathways associated with either MD resistance or susceptibility. Additionally, integrating prior information from ChIP-seq, microarray analysis, and SNPs exhibiting allele-specific expression (ASE) in response to MDV infection, we were able to provide evidence for 24 genes that are polymorphic within mEq binding sites are likely to account for gene expression in an allele-specific manner and potentially for the underlying genetic differences in MD incidence.


Introduction
Marek's disease (MD) is a lymphoproliferative disease of poultry caused by Marek's disease virus (MDV, Gallid herpesvirus 2), an oncogenic alphaherpesvirus. One of the main MD control strategies is vaccination. While MD vaccination reduces the incidence of tumor formation, it is not sterilizing, thus, does not prevent MDV from replicating or spreading amongst vaccinated birds. Additionally, field strains continue to evolve, with increased virulence in vaccinated birds. Losses with MD are further enhanced by the unpredictable and spontaneous outbreaks that occur even in vaccinated flocks [1,2]. Given the problems with vaccination, there is a need to pursue other strategies to combat MD. Identifying chickens with enhanced genetic resistance to MD is an attractive alternative to augment vaccinal control. Using genomic tools to identify genetic markers associated with MD resistance would be highly beneficial to select birds with superior disease resistance. A better understanding of the mechanisms of genetic resistance to MD would therefore contribute toward improved strategies to control the disease.
Currently, Avian Disease and Oncology Laboratory (ADOL) chicken lines 6 (MD resistant) and 7 (MD susceptible) have been developed to study the mechanisms underlying genetic resistance to MD [3][4][5]. These highly inbred (over 99%) White Leghorn lines share the same MHC haplotype, a genetic locus that has been shown to have a large effect on MD incidence [6]. Therefore, these lines enable us to focus on the remaining non-MHC genes that individually are smaller in effect size but cumulatively account for the majority of MD genetic resistance.
Selection for MD resistance is based on identifying genes where variation in their alleles is associated with variation in disease incidence. In our lab, two broad strategies have been employed and integrated. Namely, genome-wide genetic screens (e.g., QTL scans) to identify regions in the chicken genome containing the gene(s) of interest, and functional genomic screens (e.g., transcript profiling or virus-host proteinprotein interaction screens) to provide candidate genes (reviewed in 7). Despite identifying three MD resistance genes (growth hormone, SCA2, and MHC class II [8][9][10] and many more strong candidates, like other complex traits, it has been very difficult to comprehensively identify the remaining genes that are involved in resistance.
Transcriptome analysis using spleens from MD resistant and susceptible lines identified several candidate genes related to resistance and susceptibility, most of which were related to the immune response [11,12]. There are differences in the proportion of CD4 and CD8 T cells between MD resistant and susceptible lines [13], and higher expression of immunoglobulin genes in MD resistant lines when compared to susceptible lines [14]. However, these studies have not specifically examined the influence of MDV Meq, a bZIP transcription factor and the viral oncogene, and its role on genetic resistance. Using global transcriptome analysis, we identified a number of genes and pathways that are consistently associated with MD resistance or susceptibility. We also show that heterozygous SNP sites in Meq-binding sites between lines 6 and 7 are associated with allele-specific expression, which may provide a mechanism that accounts for a proportion of the variation in MD genetic resistance between these two bird lines.

Cells and culture conditions
Chicken embryo fibroblasts (CEF) from day 10 embryos were prepared from ADOL specific pathogen free line 6 and line 7, and secondary cultures plated at a density of 10 7 cells per 100 mm dish. Cells were cultured in Leibowitz's L-15 and McCoy 5A media with 15% heat inactivated fetal bovine serum, 100 U of penicillin per ml, and maintained at 37 °C. In each of two experimental replicates, three plates each that were confluent for 24 h were infected with 10 4 pfu of MDV derived from either Md5B40BAC1, our BAC clone that contains the Md5 strain MDV genome and generates virulent MDV, which we referred to as Md5 [15], or a derived nononcogenic recombinant from Md5B40BAC1 that lacks both copies of Meq, developed through recombineering [16], which we referred to as Md5∆Meq. For uninfected controls, an equal amount of uninfected CEF were added. Total RNA was isolated at 24, 48, and 96 hr from CEFs infected with either Md5 or Md5∆Meq as well as uninfected CEFs as controls.

RNA extraction microarray procedure and data analysis
Total RNA for microarray hybridization was extracted using Absolutely RNA Miniprep Kit (Stratagene, Clara, CA, USA) according to the manufacturer's instructions. RNA concentration was assessed using a Nanodrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and the RNA integrity was determined using Agilent 2100 Bioanalyzer with a RNA 6000 Nano/Pico Assay (Agilent Technologies, Palo Alto, CA). Affymetrix GeneChip Chicken Genome Arrays (Affymetrix, Santa Clara, CA, USA) were used for microarray hybridization and data collection; this chip has probe sets for annotated chicken transcripts including all 17,179 chicken unigenes and Ensembl predicted genes. RNA preparation, hybridization and scanning were performed following protocols recommended by Affymetrix by the Michigan State University Research Technology Support Facility (rtsf.msu.edu). CEL files were generated that contained the summary intensities for each probe; NCBI accession no. GSE48454. The raw data files were loaded using the Affymetrix package, and the probe intensities and normalization were done using LIMMA (linear models for microarray data) [17]. The expression value of each probe set was normalized and calibrated using the open source 'R' statistical software (version 2.11.1) through the Bioconductor project (bioconductor.org). Pair wise comparisons between the groups including Md5 vs Md5∆Meq were performed by Fishers Least Significant Difference (LSD).

Quantitative RT-PCR analysis
Total RNA was extracted using an RNeasy Mini Kit (Qiagen, Valencia, CA, USA). Complementary DNA was synthesized using a Superscript III First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA). Gene expression levels were measured using SYBR Green PCR Master Mix (Invitrogen) on an ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The β-actin gene was used for normalization and each target gene was analyzed in triplicate. The PCR conditions were as follows: 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 sec, 60 °C for 1 min, and 95 °C for 15 sec. At the end of amplification, a melting curve analysis was done by heating the PCR products to 65-95 °C, which was held for 15 sec at increments of 0.2 °C to measure fluorescence and confirm the presence of a single amplification product. For negative controls, no-RT was used as template in place of single-stranded cDNA in the qRT-PCR. The data analysis was performed with the comparative ∆∆Ct relative quantification method [18].

Pathway analysis
The differentially expressed genes (DEGs) were identified amongst uninfected, Md5-infected, and Md5∆Meq-infected lines 6 or 7 CEF. These genes were further analyzed for inclusion in Gene Ontology (GO) categories and pathways in order to examine their biological processes. Categorization of genes based on significant biological properties was done using Gene Ontology Project (http://www.geneontology.org/). The genes were grouped categories based on common biological properties.
The pathway analysis was carried out using Ingenuity Pathway Analysis software (IPA, Ingenuity Systems, Redwood City, CA, USA). The annotated genes were grouped into networks, functions, or canonical pathways. The data with gene IDs and expression fold-change were uploaded into the software. The gene IDs were mapped into its corresponding gene object in the Ingenuity Knowledge Base (IKB). The network of focus genes were generated based on the information contained in the IKB into a global molecular network. Functional analysis generates biological functions that are significant to the genes in the data uploaded. The canonical pathways are generated based on the DEGs from the data. All the IPA analyses were carried out based on Fisher's exact test to determine the association between the differentiallyregulated genes and the network, biological function, and canonical pathway.

Identification of Meq-dependent genes related to genetic resistance and susceptibility
To determine the role of Meq in MD genetic resistance, we performed global transcriptome analysis in CEFs from lines 6 and 7, two highly inbred experimental White Leghorn chicken lines that differ greatly in MD incidence. To determine if genes in Meq-regulated pathways are influenced by the genetic resistance status of the host, the CEFs from both chicken lines were infected with either Md5, which is a fully virulent MDV, or Md5∆Meq, a recombinant MDV that lacks both copies of Meq and is nonocogenic. We chose three time points: 24, 48 and 96 hr to study the gene expression changes induced during these phases of virus infection. We hypothesized genes that are uniquely differentially expressed in Md5-infected groups and not in Md5∆Meq-infected cells indicate ones that are directly or indirectly dependent on Meq expression.
We made pair-wise comparisons to generate a list of genes that are unique to lines 6 and 7, which were further classified into genes dependent or not dependent on the expression of Meq as shown in (Figure 1). The details on DEGs at each time point in Md5-infected and Md5∆Meq-infected groups for each line are provided in Table 1. GO categorization for pathway analysis of lines 6 and 7 were performed using a 'union' of DEGs in all three time points. The number of Meq-dependent genes that are associated with MD resistance and susceptibility are shown in Figure 2A and 2B, respectively. The DEGs detected by microarrays were validated by qRT-PCR for a randomly subset of genes that were selected from the top biological networks. The fold-change as determined by qRT-PCR was significantly correlated to the findings from microarray data (r 2 =0.67; P<0.05) ( Figure 3).
To further interpret if genes in Meq-regulated pathways are influenced by the genetic resistance status of the host, we combined information about Meq-bound and regulated [19] with the list of DEGs from lines 6 and 7 at a significance level of P<0.01. Interestingly, about 20% of the genes from line 6 and 35% of the genes from line 7 overlapped with Meq-bound and regulated genes previously identified in vitro ( Figure 4).

GO categorization of Meq-dependent gene list in CEF from MD resistant or susceptible chickens
The DEGs of Md5-infected and Md5∆Meq-infected groups in each of lines 6 and 7 were classified into different functional GO categories. Each of the putative genes associated with MD resistance or susceptibility were assigned to molecular function categories as designated by the GO database. Based on functional annotation clustering using the highest classification stringency, there were 53 and 44 clusters in DEGs related to MD resistance and susceptibility, respectively.  The most significant GO categories of differentiallyexpressed Meq-dependent genes associated with MD resistance with enrichment cutoff set to P<0.05 ( Figure 5A) were transcription regulator activity, inflammatory cell apoptosis, immune response, positive regulation of apoptosis, and negative regulation of cell proliferation whereas differentially-expressed Meq-dependent genes associated with MD susceptibility ( Figure 5B) were sequence-specific DNA binding, regulation of cell proliferation, DNA-dependent regulation of transcription, and protein heterodimerization. The significant categories that were common to both line 6 and line 7 Meq-dependent genes were transcription activator activity and transporter activity.
To further examine the biological pathways that are altered during the host response to MDV infection in the presence of Meq, the DEGs associated with MD resistance and susceptibility were analyzed. The DEGs associated with MD resistance that were Meq-dependent were significantly associated (P<0.05) with 21 canonical pathways ( Figure 6A). Some of the top pathways included apoptosis, death receptor signaling pathway, and Myc-mediated apoptosis, which plays an important role for controlled cell death. Genes involved in cell death were up-regulated and genes involved in cell maintenance were down-regulated in this category. The pathways involved in DEGs associated with MD resistance that were not dependent on Meq included DNA replication pathway,     Similar pathway analyses were performed on the DEGs associated with MD susceptibility that were Meq-dependent (19 significant canonical pathways; Figure 6B) and Meqindependent groups. Some of the top pathways with DEGs associated in MD susceptibility that were Meq-dependent group included cell cycle regulation, MAPK signaling pathway, and Jak-STAT signaling pathway. The pathways involved in DEGs associated in MD susceptibility that were not dependent on Meq included VEGF signaling pathway, DNA replication pathway, Hedgehog signaling pathway, mismatch repair, and insulin signaling pathway. Of these analyses, there were three pathways that were common between the two groups.

Allele-specific binding bias of Meq at SNPs between line 6 and line 7
We have previously identified SNPs exhibiting allele-specific expression (ASE) in response to MDV infection using F 1 birds derived by intermating lines 6 and 7 [20]. However, the mechanism for the allelic transcriptional variation in response to MDV infection is unknown. Thus, we explored whether there might be association of SNPs in the promoters of specific genes with variation in Meq-dependent binding. Specifically, we screened for genes that (1) are directly regulated by Meq, (2) exhibit ASE in response to MDV infection, and (3) have a SNP between lines 6 and 7 in the Meq-binding motif. A total of 24 candidate genes exhibited polymorphisms between the two lines in their Meq-binding motifs (e.g., TGACTCA or CACACAGC) located in their proximal promoter region (Table  2). Interestingly, this analysis revealed genes in the MAPK kinase signaling pathway (e.g., BRAF and GPS1), contributing to cell cycle regulation (e.g., FANCA, NDC80, and RB1), or coding for transcription factor subunits (e.g., AP1M1 and OPTN), all of which are known to be connected to response to virus infection and members of pathways previously implicated in our prior studies [19,20].

Discussion
In recent years, unpredictable and spontaneous vaccine breaks resulting in devastating losses to poultry farms [1], have further necessitated the need to explore alternate strategies for MD prevention. Selection for increased genetic resistance to MD is a control strategy that has proven to work and can augment MD vaccine protection [21]. Over the years, several attempts have been made to identify candidate genes that determine genetic resistance to MD. Based on a variety of genetic and genomic strategies, these studies have shown various factors underlying the mechanism of resistance and susceptibility [5,8,12,[21][22][23][24]. It is imperative to further understand the underlying mechanisms of MDV pathogenesis and further examine how this impacts genetic resistance. Being a complex disease, variability in a single gene cannot explain the basis for genetic resistance by itself. Although some previous studies have examined differential gene expression patterns after exposure to MDV in MD resistant and susceptible lines of chicken, there are no studies that have explored the role of Meq in determining genetic basis for resistance. Among several viral genes, null mutants for Meq showed no oncogenicity while knock out mutants of other viral genes like pp38 and ICP4 only resulted in attenuated virulence [25. 26]. Given the importance of Meq, we attempted to provide insights into the molecular mechanisms of MD resistance in the presence of Meq in response to MDV infection.
Through selective breeding and identification of phenotypic variation with respect to MD incidence, genetically resistant and susceptible lines of chicken have been developed [5]. One of the major contributors to resistance is variable MHC haplotypes, as evidenced by differential MHC haplotypes responsible for phenotypic variation in lines N and P [27]. However, there are numerous other factors that are non-MHC dependent which play an important role in genetic resistance to MD. Lines 6 and 7 share the same MHC haplotype but vary markedly in their susceptibilities to MD [3,28]. Hence we used this model to unravel non-MHC related basis for variable susceptibility to MD. Moreover, our study adds an essential aspect that has not previously been explored by using MDV constructs with and without Meq, its principal oncogene.
Further comparisons of results from our previous findings on genome-wide regulatory network of Meq [19] should provide additional confidence in the declared list of genes influencing resistance and susceptibility. This helps us to determine if genes in Meq-regulated pathways are influenced by the genetic resistance status of the host. We found that more than 20% of genes from line 6 and 35% of genes from line 7 were overlapping with genes that had binding sites for Meq and were transcriptionally regulated by Meq.
We hypothesized that the Meq-dependent differentially expressed genes would be involved in the downstream molecular pathways that might play an important for maintaining MD resistance. In order to identify the types of specific molecular functions of these genes, the Meqdependent differentially expressed genes were annotated using the GO. Each gene was assigned to Molecular Function and Biological Process categories as designated by the GO database. The molecular functions and the biological process were significantly different between MD resistant and susceptible lines. The resistant line was enriched for positive regulation of cell death whereas the susceptible cell line was enriched for regulation of cell proliferation. Further analyzing the cellular pathways involved, we found that apoptosis signaling and death receptor pathways were among the significantly represented pathways. Specifically, we found an up-regulation of caspases 3, 6, and 8. Caspase 3 is an executioner caspase that is activated by both extrinsic and mitochondrial intrinsic pathways of apoptosis. Caspase 3 is activated by caspase 8 and activated caspase 3 can in turn activate other caspases like 6 and 7, which eventually leads to cell death [29]. We also noted a down-regulation of inhibitors of apoptosis like BIRC2 and antiapoptotic proteins like Bcl-2 and Bcl-X L . Interestingly, we have shown in a previous study that most of these genes have Meq-binding sites and are transcriptionally regulated by Meq [19]. The up-regulation of these caspases and down-regulation of anti-apoptotic factors could be one of the determinants of MD resistance in line 6. We found that mitogenic signals like MAPK signaling and regulation of cell cycle involving cyclin D were overrepresented in the line 7 transcriptome. We have previously shown that one of the major cellular pathways that Meq targets to induce transformation is the MAPK signaling pathway. Specifically, we have shown that Meq up-regulates mitogenic signals like MEK1, MEK2 and Ras, which drive the cells towards proliferation. Also, Meq down-regulates inhibitory signals like phosphatases that limit the activation of MAPK signaling. In the present study, we found that at least some of these genes are associated in differential MDV resistance/susceptibility. In the MD resistant chicken line, there was a down-regulation of MAP2K2 (MEK2) in Md5-infected group compared to Md5ΔMeq-infected group. This highlights the role of Meq in modulating MAPK signaling and how down-regulation of mitogenic signals may be associated in genetic resistance to MD. In contrast, we found an up-regulation of Ras, another mitogenic signal in line 7 infected with Md5 compared to the Md5ΔMeq-infected group, raising the possibility that upregulation of mitogenic signals underlies genetic susceptibility to MD. These results suggest that the genes related to cell death and proliferation is one of the major determinants of genetic variability in resistance/susceptibility to MD. In one of our previous studies using RNA sequencing analysis, we have shown that apoptotic signaling is significantly overrepresented in spleens of birds infected with MDV [30]. Our group has previously demonstrated that Jak-STAT is an important cellular pathway involved in genetic basis of MD resistance [20,30]. Similarly, we also noted that genes in Jak-STAT pathway are transcriptionally up-regulated in susceptible line 7 and not in the resistant line 6. Activation of Jak-STAT pathway results in nuclear translocation of activated STAT dimer resulting in transcription of genes involved in cell survival and proliferation [31]. In addition to MAPK pathway, this could be another mitogenic signaling pathway that is involved in transformation by Meq.
We have previously found a subset of transcription factor binding sites from ChIP-seq analysis. Further comparison with our previous findings on the ASE between lines 6 and 7 provides additional confidence in the list of genes regulated by Meq that exhibit polymorphisms at the proximal promoter region and allele-specific gene expression. Additionally, we also found that many allele-specific binding shows a bias between line 6 and line 7, which potentially represents the functional distinction between MD resistance and susceptibility. If planned experiments show an association between these SNPs and MD incidence, then it would validate the importance of Meq-regulated transcription in MD genetic resistance.
We also found down-regulation of IRG1, a putative proapoptotic factor which was recently described as a candidate susceptibility gene for MD [11] in line 7 but no expression in line 6. The fact that we noted this transcriptional response in line 7 infected with Md5 and not Md5∆Meqinfected group raises an important and novel corollary that Meq has a role in regulating IRG1 expression. Further, we also noted transcriptional regulation of other genes (STAT1, MyD88, IFN-γ) proposed in the biological interaction network analysis of IRG1 in the previous study. In addition to corroborating the results from that study, this further underscores the role of Meq in modulating apoptotic signaling as a determinant of MD susceptibility.
Admittedly, CEFs are not the natural target cells for MDV. However, it is worth noting transcript profiling using spleens and other organs suffers from the potential confounding influence of differences in cellular composition and reagents to properly quantify each cell subtype are limiting in chicken. Thus, the high degree of overlap between our current results, which utilized cultured cells from defined chicken lines, and those from our other studies that used actual birds [11,12,20,30] for the identification of critical pathways suggests that CEF are reasonable models. This agreement also suggests an important contribution of innate immunity towards MD genetic resistance. Furthermore, as cells and viruses can be more easily manipulated and monitored compared to actual MDV challenges of live birds, this suggests that more precise experiments or hypotheses should be first explored in vitro prior to animal studies, especially when one wishes to screen for the direct influences of particular viral proteins like Meq.
In conclusion, we have made significant insights on the different sets of genes and pathways that interact to modulate MD resistance and susceptibility. Taken together, our findings add to the current understanding of the mechanism behind Meq-induced responses that lead to MD resistance or susceptibility. The overlap between SNPs from ASE and ChIPseq data show that allele-specific Meq binding at many of these SNPs could allow a functional distinction between line 6 and line 7. We have confirmed that allele-specific binding biases by Meq tend to be one of the underlying genetic mechanisms for line 6 and line 7 alleles. In addition, this study forms the basis for selection of candidate genes that might be associated in genetic resistance to MD.