RNA Sequencing revealed differentially expressed genes functionally associated with immunity and tumor suppression during latent phase infection of a vv + MDV in chickens

Very virulent plus Marek’s disease (MD) virus (vv + MDV) induces tumors in relatively resistant lines of chickens and early mortality in highly susceptible lines of chickens. The vv + MDV also triggers a series of cellular responses in both types of chickens. We challenged birds sampled from a highly inbred chicken line (line 63) that is relatively resistant to MD and from another inbred line (line 72) that is highly susceptible to MD with a vv + MDV. RNA-sequencing analysis was performed with samples extracted from spleen tissues taken at 10-day and 21-day post infection (dpi). A total of 64 and 106 differentially expressed genes was identified in response to the vv + MDV challenge at latent phase in the resistant and susceptible lines of chickens, respectively. Direct comparisons between samples of the two lines identified 90 and 126 differentially expressed genes for control and MDV challenged groups, respectively. The differentially expressed gene profiles illustrated that intensive defense responses were significantly induced by vv + MDV at 10 dpi and 21 dpi but with slight changes in the resistant line. In contrast, vv + MDV induced a measurable suppression of gene expression associated with host defense at 10 dpi but followed by an apparent activation of the defense response at 21 dpi in the susceptible line of chickens. The observed difference in gene expression between the two genetic lines of chickens in response to MDV challenge during the latent phase provided a piece of indirect evidence that time points for MDV reactivation differ between the genetic lines of chickens with different levels of genetic resistance to MD. Early MDV reactivation might be necessary and potent to host defense system readiness for damage control of tumorigenesis and disease progression, which consequently results in measurable differences in phenotypic characteristics including early mortality (8 to 20 dpi) and tumor incidence between the resistant and susceptible lines of chickens. Combining differential gene expression patterns with reported GO function terms and quantitative trait loci, a total of 27 top genes was selected as highly promising candidate genes for genetic resistance to MD. These genes are functionally involved with virus process (F13A1 and HSP90AB1), immunity (ABCB1LB, RGS5, C10ORF58, OSF-2, MMP7, CXCL12, GAL1, GAL2, GAL7, HVCN1, PDE4D, IL4I1, PARP9, EOMES, MPEG1, PDK4, CCLI10, K60 and FST), and tumor suppression (ADAMTS2, LXN, ARRDC3, WNT7A, CLDN1 and HPGD). It is anticipated that these findings will facilitate advancement in the fundamental understanding on mechanisms of genetic resistance to MD. In addition, such advancement may also provide insights on tumor virus-induced tumorigenesis in general and help the research community recognize MD study may serve as a good model for oncology study involving tumor viruses.

www.nature.com/scientificreports www.nature.com/scientificreports/ Differentially expressed genes in response to MDV challenge. To improve the reliability and comparability of differential expression analysis, a total of 10,257 identified genes with FPKM values ≥1 in all treatment groups were included in this and subsequent analyses. At 10 dpi, 61 and 79 DEGs in response to MDV challenge were identified in the line 6 3 and 7 2 birds, respectively. Of which, 54 (88.52%) and 46 (58.22%) DEGs were significantly up-regulated in the line 6 3 and 7 2 birds, respectively (Fig. S1A,C). At 21 dpi, 11 and 60 DEGs in response to MDV challenge were identified in the line 6 3 and 7 2 birds, respectively. Interestingly, all the 11 DEGs of line 6 3 were upregulated while 51 (85%) of the 60 DEGs in line 7 2 were upregulated ( Fig. 2A,B, and Fig. S1B,D). Detailed lists of all DEGs induced by MDV challenge in the two genetic lines at 10 and 21 dpi are given in Supplementary Tables S2-S5. As shown in the Venn diagrams (Fig. 2C,D), there were dozens of mutually exclusive DEGs between the lines 6 3 and 7 2 at 10 and 21 dpi except the line 6 3 at 21 dpi with only a single mutually exclusive DEG, in addition to 26 and 10 DEGs identified in common between the lines of birds at 10 and 21 dpi, respectively. Within each of the lines between the two time points, 10 and 21 dpi, there were a few (2 in line 6 3 at 21 dpi) up to dozens of DEGs identified mutually exclusive between the two time points, in addition to 9 and 33 DEGs in common between the two time points in line 6 3 (19) and mutually exclusive (55 and 29) between the line 6 3 and 7 2 control groups' contrast and the MDV challenged groups' contrast, respectively at 10 dpi. (D) A Venn diagram showing the number of DEGs in common (18) and mutually exclusive (24 and 78) between the line 6 3 and 7 2 control groups' contrast and the MDV challenged groups' contrast, respectively, at 21 dpi. (E) A Venn diagram showing the numbers of DEGs in common and mutually exclusive between the line 6 3 challenged/control, line 7 2 challenge/control, and the challenged groups of lines 6 3 and 7 2 at 10 dpi. (F) A Venn diagram showing the numbers of DEGs in common and mutually exclusive between the line 6 3 challenged/control, line 7 2 challenge/control, and the challenged contrast to line 7 2 control groups, 36 and 32 DEGs were significantly expressed at higher levels, 38 and 10 DEGs were expressed at lower levels in the line 6 3 control groups at 10 and 21 dpi, respectively (Fig. 3A, Fig. S2A,B, and Tables S6, S7). Following MDV challenge, 32 and 59 DEGs were expressed at significantly higher levels, 16 and 37 DEGs were expressed significantly at lower levels in the line 6 3 compared to the line 7 2 at 10 and 21 dpi, respectively (Fig. 3B, Supplementary Fig. S2C,D, Supplementary Tables S8, S9). There were dozens of genes expressed differentially between the two lines regardless of MDV-challenge, and there were 19 and 18 genes differentially expressed in common between the two lines as well at 10 (Fig. 3C) and 21 dpi (Fig. 3D), respectively. In MDV-challenged birds, there were certain degrees of overlaps of genes differentially expressed between the two lines at 10 (Fig. 3E) and 21 dpi (Fig. 3F). A close examination of 5 and 12 genes at 10 and 21 dpi, respectively, between the lines 6 3 and 7 2 for both the control groups and the MDV challenged groups is depicted in Fig. 3G,H. The expression levels of those genes were significantly differentiated between the two lines, but none of the expression levels of those genes was between the control groups. ddPCR validation of the RNA-Seq data. To assess the validity of RNA-Seq data, a total of 12 pairs of primers (Table S10) was designed targeting 12 selected DEGs that were identified in a combination of 23 comparisons of this study by time point (dpi), MDV treatment and chicken lines. These primers were used in droplet digital TM PCR (QX200TM ddPCR system; Bio-Rad Laboratories, Inc., Hercules, CA, USA) analysis for the validation. Three biological replicates were used for each treatment group. A scatterplot was generated by comparing the log 2 FC determined by RNA-Seq analysis and ddPCR analysis. The result showed there was a fairly high correlation between the two groups of log 2 FCs determined by the two methods (Fig. S3, R 2 = 0.496, P value < 0.001). The ddPCR data validated the RNA-Seq data. ddPCR analysis of MDV microRNAs' expressions associated with MDV latency. The expressions of two MDV microRNAs, miR-M3-3p and miR-M12-3p, in each of the total RNA samples extracted from MDV challenged birds and included in the current RNA_Seq analysis were quantitatively assessed by ddPCR to ensure that each of those total RNA sample of this study was from a challenged bird that had entered the latent phase. The miR-M3-3p and miR-M12-3p have been characterized with significantly higher expression in latent phases in contrast to that at the early cytolytic phase 39 . The ddPCR analysis was conducted with technical replicates using customer primers (Table S11). Splenic total RNA samples extracted from MDV challenged birds at 5 DPI were also used to contrast the difference of the microRNA expressions. The 5-DPI total RNA samples were aliquoted from a sister project conducted simultaneously along with this project using the same hatch of line 6 3 and line 7 2 birds under the same exact conditions. In a 20 μL reaction, an average of 7,020 ± 653.9 and 4,587 ± 610.0 miR-M3-3p copies at 10 DPI and 5,123 ± 866.6 and 3,750 ± 1,092.7 copies at 21 DPI was detected for lines 6 3 and 7 2 , respectively; in contrast, an average of 422 ± 30.6 and 356 ± 35.9 copies was detected in the 5 DPI samples of line 6 3 and 7 3 , respectively. The average of miR-M3-3p expression copies among all the 10 and 21 DPI individual total RNA samples ranged from 2,760 ± 180 to 8,120 ± 40 for line 6 3 , and 1,860 ± 20 to 6,470 ± 470 for line 7 2 .
The results of Analysis of Variance showed both MDV microRNAs' expressions at 10 and 21 DPI significantly differed from those at 5 DPI (P < 0.01; leverage plots of the MDV microRNA expression data are given in Figs. S4 and S5 graphically illustrating the differences in expression at 5, 10, and 21 DPI). This in combination with the ranges of individual bird MDV microRNA expression evidently suggested that none of the MDV challenged birds under this study was remained at the early cytolytic phase.

Functional analysis of the DEGs.
To better understand the MDV-induced DEGs, functional enrichment analysis was conducted for four separate sets of the DEGs, which included two sets of MDV-induced DEGs in lines 6 3 and 7 2 at 10 dpi and two sets of MDV-induced DEGs in lines 6 3 and 7 2 at 21 dpi. At 10 dpi, the up-regulated genes in response to MDV challenge in line 6 3 were primarily enriched in 13 GO terms, which are primarily associated with the immune system including immune effector process, immune system process, immune response, defense response, and regulation of immune response. The down-regulated genes in line 6 3 at 10 dpi were associated with cell communication (FAM132 A, FN1 and PROKR2) and immune system process (AQP3). The up-regulated genes of line 7 2 at 10 dpi were also significantly enriched in two GO terms of immune response and immune effector process. Some up-regulated DEGs of both lines 6 3 and 7 2 were also involved with the Influenza A pathway. In contrast, the line 7 2 down-regulated genes were enriched in four GO terms and four KEGG pathways. Notably, these functional categories included several terms involved in fatty acid metabolism, such as fatty acid groups of lines 6 3 and 7 2 at 21 dpi. (G) Depicting both the expression levels of 5 genes at 10 dpi and (H) 12 genes at 21 dpi without (control) and with MDV challenge between the line 6 3 and line 7 2 groups. A red star indicates a gene that was differentially expressed (FDR <0.05 and FC >2) in response to MDV challenge in the line 7 2 birds, whereas a green star indicates that gene was differentially expressed (FDR <0.05 and FC >2) in response to MDV challenge in the line 6 3 birds. Both groups of charts in G and H demonstrated that the difference in gene expression in the control groups did not alter the differential expression status of genes in response to MDV challenge, at least in this study. (2019) 9:14182 | https://doi.org/10.1038/s41598-019-50561-x www.nature.com/scientificreports www.nature.com/scientificreports/ biosynthetic process, biosynthesis of unsaturated fatty acids and fatty acid metabolism (Table 1). Of the 11 significantly up-regulated genes in line 6 3 at 21 dpi (Table S4), two were granzyme (GZMA and GZMK) genes, three were immune response genes (MX, OAS*A and CCLI10), one was an Interleukin Four Induced gene 1 (IL4l1) and one was avidin gene (AVD). For line 7 2 , the up-regulated DEGs were significantly enriched in GO terms primarily associated with immune response, including defense response, innate immune response, and response to interferon-gamma, in addition to influenza A pathway ( Table 2).
Gene function analysis showed that DEGs that were relatively expressed at higher levels in the unchallenged line 6 3 birds than those of the unchallenged line 7 2 birds were significantly enriched in extracellular structure organization and extracellular matrix organization, as well as pathways of Focal adhesion and ECM-receptor interaction ( Table 3). DEGs that were relatively expressed higher in the unchallenged line 7 2 birds than those of the unchallenged line 6 3 birds were mainly involved in metabolic pathways, such as pathways of Butanoate metabolism, Fatty acid metabolism, Biosynthesis of amino acids. Following infection at 10 dpi, six of the DEGs were immune genes (GAL1, GAL2, GAL7, ABCB1, LEPR and SNED1) and two were involved in NOD-like receptor signaling pathway (K60 and HSP90B1), all were upregulated in expression in line 6 3 birds. In contrast, four DEGs associated with response to stimulus (HSP90AB1, GNG4, ANXA1, CD180) were upregulated in expression in line 7 2 birds. At 21 dpi, upregulated DEGs identified in line 6 3 were significantly enriched in GO terms involved in system development and two pathways including Focal adhesion and ECM-receptor interaction, while upregulated DEGs identified in line 7 2 were mainly associated in immunity, as indicated in enriched functional categories such as antigen processing and presentation, adaptive immune response, cell killing and KEGG pathway (Table 4).

Key genes for MD resistance.
To further identify key genes that may confer genetic resistance to MD from all the unique DEGs, we focused on the followings: (1) DEGs located within reported QTL regions of MD resistance (http://www.animalgenome.org/cgi-bin/QTLdb/GG/index); (2) DEGs with large fold change (FC >2) in one line but not in the other (FC < 0.5); (3) DEGs differentially expressed between the two lines regardless of MDV challenge; (4) MDV induced DEGs with differential expression between infected line 6 3 and line 7 2 birds (Fig. 3E-H). These comparisons resulted in a total of 104 unique DEGs. Each of these DEGs was subjected to a separate Gene Ontology (GO) analysis and a comprehensive literature review. Jointly considering the expression pattern of each gene in different chicken lines and the potential functionality, a short list of 27 top DEGs was selected as the likely candidate genes identified from this study, which may potentially play key roles in conferring genetic resistance to MD (Fig. 4). These genes were broadly grouped into three functional groups, including virus process, immunity and tumor suppression.

Discussion
Several studies have documented the global gene expression responses to MDV infection in chicken 3,12,28-36 . These analyses, while providing valuable and important mechanistic clues to host responses and virus-host interactions, are of limited scope due to coverage of gene arrays used and lack information of comparisons between chickens exhibiting different resistance to MD. By this study, we further expanded the current knowledge with regard to host transcriptomic response upon MDV challenge through comparing the whole transcriptome changes between a highly inbred MD-resistant and -susceptible lines of chicken using RNA-Seq technology. This study was primarily focused on the latent stages post MDV inoculation. Our results revealed quite different transcriptome patterns, which may provide new insights potentially conferring genetic resistance to MD.
Close examinations of the transcriptomic changes in response to MDV challenge showed that many genes associated with immunity and anti-tumor functions were activated in both lines of chickens. For example, two granzyme genes, GZMA and GZMK, that are associated with apoptosis, and the MX1 gene that is a well-known IFNs-induced gene 40 , were up-regulated in both lines at 10 and 21 dpi, suggesting on going interactions that took place between host and the MDV pathogen. In contrast, more different features in response to MDV infection between the two inbred lines were observed. More DEGs were identified in the MD susceptible line 7 2 than those in the resistant line 6 3 at both 10 and 21 dpi ( Fig. 2A), indicating that MDV induced much stronger response at the transcriptional level in line 7 2 than in line 6 3 birds. However, it is interesting to note that the proportion of up-regulated DEGs was much higher in line 6 3 than that in line 7 2 (Fig. 2B), which suggested that insufficient immune response was probably activated against tumorigenesis regardless more immune-related gene expressions were altered in response to MDV challenge. Specifically, in MD resistant line 6 3 , a notable immune response was observed at 10 dpi, which was characterized by the overrepresentation of DEGs enriched in many immune response-related GO terms, such as immune response, defense response to virus, positive regulation of immune response, and so on. While at 21 dpi, only 11 genes were identified to exhibit a statistically significant change in expression, indicative of a minimal change compared to normal tissues. In susceptible line 7 2 , a totally different transcription pattern was revealed across the two-time points. At 10 dpi, 33 out of 79 DEGs (41.77%) were down-regulated in expression (Fig. 2B, Fig. S2C). These down-regulated genes were significantly enriched in multiple GO terms and pathways, including several terms associated with fatty acid metabolism, such as fatty acid biosynthetic process, pathways of fatty acid metabolism, and biosynthesis of unsaturated fatty acids (Table 2). Genes involved in these terms included SCD, FADS1, FADS2 and HPGDS. Previous studies showed that Hepatitis C virus (HCV) could induce imbalance in lipid homeostasis in host cells 41,42 . However, it seemed that inhibition of fatty acid biosynthesis can suppress virus replication 43,44 and control cancer cell proliferation 45 . Therefore, whether the down-regulation of these genes in susceptible line 7 2 playing a role in establishment of latency or representing a host response that may contribute to repressing the activities of MDV infection needs to be further investigated. At 21 dpi, 60 DEGs were identified being differentially expressed compared to the control birds and 51 of those were up-regulated. Gene function analysis showed that the up-regulated genes were significantly enriched in immune response-related terms such as innate immune response, interferon-gamma-mediated signaling pathway, and so on. In particular, the IRF1 (Interferon regulatory factor-1) gene expression was upregulated by 4.2 folds (from an FPKM of 1278.1 to 5415.7). IRF1 is a member of the IRF gene family of transcription factors that  www.nature.com/scientificreports www.nature.com/scientificreports/ binds to the virus-inducible cis-elements of IFN-α and IFN-β gene promoter as well as to the interferon response sequence of IFN-inducible gene promoters 46 . Accordingly, the 51 identified up-regulated DEGs were significantly enriched in IRF-1 (P value = 5.66E-5) by a motif enrichment analysis using the g:Profilier (http://biit.cs.ut. ee/gprofiler/). A total of 15 up-regulated genes containing potential IRF-1 binding site was revealed, including IFIH1, BCL2L14, GZMK, ISG12(2), SERPINB10, IFI27L2, AVD, TAP1, CCL19, CMPK2, ENSGALG00000026152, ENSGALG00000001629, ENSGALG00000006384, ENSGALG00000013057 and ENSGALG00000019141. These gene functions implied strong defense responses might have been activated in susceptible line 7 2 at 21 dpi, which is in good agreement with earlier observations that MDV were reactivated in susceptible line of birds at late latency stage 1 . Furthermore, we observed some genes showed extreme difference in response to MDV challenge between the two genetic lines of birds at each time points (Fig. S2C,D). These genes might also play key roles modulating genetic resistance to MD, therefore, follow-up investigations are warranted. One of those DEGs is IL4I1. Reportedly it is an inhibitor of the CD8+ antitumor T-cell response and may facilitate tumor growth. Our data showed that the expression of IL4I1 was significantly increased in the susceptible line 7 2 birds at 10 dpi, while it remained unchanged in the resistant line 6 3 birds at the same time point, which indicated that genes like IL4I1 may contribute to immunosuppression and facilitate tumorigenesis through inhibiting the CD8+ antitumor T-cell response in MD susceptible chickens. If so, it is then in good agreement with the functional findings of this gene in mice 47 .
Differences in gene expression levels were also observed between the two lines at both time points for birds without and with MDV challenge. Notably, genes associated with Focal adhesion and ECM-receptor interaction were observed with higher expression in the resistant line 6 3 control birds (Table 3). Interestingly, following MDV challenge, genes involved with these two pathways exhibited consistently elevated expression in the resistant line 6 3 birds than that in the susceptible line 7 2 birds at 21 dpi (Table 3), indicative of potential roles of these two pathways in MD resistance. The Focal adhesion and ECM-receptor interaction pathways have a profound influence on major cell programs including growth, differentiation, migration, and survival 48 , which suggests differences in functional properties of cells between the resistant and the susceptible lines of chicken. Furthermore, earlier studies have documented that these two pathways play a major role in immunity 49,50 . More interestingly, these pathways are also reportedly linked to tumor progression 51,52 . Therefore, the higher expression of these genes involved in Focal adhesion and ECM-receptor interaction may be contributing to tumor formation and development. Notably, the majority of the significantly enriched functional categories for DEGs with higher levels of transcriptional expression were observed in the susceptible line 7 2 birds other than in the resistant line 6 3 birds in response to MDV challenge, which are related to immune response, including positive regulation of T cell mediated immunity, adaptive immune response, and positive regulation of cell killing. This, again, implicated the strong reactivation of the susceptible line 7 2 birds to MDV infection at these points of time.
Merging the DEGs identified between the infected and uninfected birds as well as between the two lines resulted to a total of 284 unique genes. Among them, 189 DEGs (66.6%) have been previously reported in studies of host response to MDV challenge 3,12,[28][29][30][31][32][33][34][35][36] , including genes associated with chemokine (CCL1, CXCL13L2, CXCL12 and CCL19), cytokine (IL21R, IL2RA, LEPR, IL12RB2, FLT3), innate immune response (IRF1, STAT1,  CD36, IFIH1, SLC11A1, IGJ, PAPR9, GCH1, LAG3, RSAD2 and HPX), and adaptive immune response (B2M,  P2RX7, SLC11A1, ANXA1, RSAD2, HPX, and BFIV21). This overlap confirms the roles of these candidate genes  www.nature.com/scientificreports www.nature.com/scientificreports/ that may play important roles in response to MDV challenge. On the other hand, due to other differences including genetic background of chickens, MD virus strains, tissue samples and sampling time, novel DEGs were identified in this study, which expanded the candidate gene pool conferring genetic resistance to MD. Jointly considering the gene expression pattern in two lines, gene function, and comparison with QTL regions, 27 promising top candidate genes were proposed based on this study, which may highly likely play key roles in MD including conferring genetic resistance to MD (Fig. 4). These genes showed large differences either between resistant and susceptible lines, or in response to MDV challenge within each of the genetic lines, or both, or are located at reported MD-related QTL regions. Functional analysis classified these genes into three broadly functional categories, including viral process, immunity, and tumor development.
Earlier studies have shown that viral load was different between the line 6 3 and line 7 2 at latency phase of infection 25,26 . Consistent with these observations, we found two genes, F13A1 and HSP90AB1, with known functions   53,54 , that were significantly induced in expression in response to MDV challenge in the susceptible line 7 2 birds, but with little change in the resistant line 6 3 birds at 10 dpi. Presumably, these two genes might be hijacked by MDV to facilitate viral replication in the susceptible line 7 2 birds. Additionally, there were 19 genes that may potentially play roles in immunity. Among them, ABCB1LB, RGS5, C10ORF58, OSF-2, MMP7 and CXCL12 were observed with higher expression levels in the resistant line 6 3 birds regardless of MDV challenge treatment. Our results also revealed that some immune response-related genes, including HVCN1 and PDE4D, were down-regulated in the susceptible line 7 2 birds in response to MDV challenge. Together with the identification of up-regulation of IL4I1 and down-regulation of defensing genes, GAL1, GAL2 and GAL3, in the susceptible line 7 2 , it is postulated here that the susceptible line 7 2 birds may suffer from immunosuppression at latent stage of MDV infection, which is in good agreement with a previous report 29 . Other promising candidates involved in immunity were uniquely induced by MDV in the resistant line 6 3 (PARP9, EOMES, MPEG1, PDK4, CCLI10, K60 and FST) compared to the susceptible line 7 2 . Presumably, these genes may belong to the group that contribute to genetic resistance to MD. MD is a virus-induced tumorous disease of chicken and has been proposed to be an invaluable model for investigation of virus-induced cancers 7,55,56 . Our results showed that six putative tumor suppressor genes, including ADAMTS2 57 , LXN 58 , ARRDC3 59 , WNT7A 60 , CLDN1 61 and HPGD 62-64 , exhibited higher expression levels in the resistant line 6 3 than in the susceptible line 7 2 birds, which should confer genetic resistance to MD. The findings from this study underscored the value of MD study serving as a model for better understanding tumorigenesis and raised the possibility that more novel genes involved in tumorigenesis and development need to be further explored.

Conclusion
In summary, we have carried out a comprehensive gene expression study using RNA-Seq analysis and identified hundreds of differentially expressed genes in response to a vv + MDV challenge in two highly inbred lines of chickens at 10 and 21 days post MDV infection. It is well-documented that one of the lines, the line 6 3 , is relatively resistant to MD, while the other, the line 7 2 , is highly susceptible to MD. We identified a total of 284 unique coding genes that likely affect the resistance to MD, which provided a sizable and valuable addition to the current candidate gene pool that is reportedly involved with genetic resistance to MD. We showed that the response pattern of gene expression at 21 dpi supported a reactivation of MDV in the susceptible line. We further proposed 27 promising candidate genes that may play key roles conferring genetic resistance to MD. Notably, most of these promising candidate genes identified in this study are reportedly associated with immunity and tumor suppression, which may have an important implication on virus-induced tumorigenesis in general and highlighted the value of MD model for tumor virus-induced cancer study. Further directions should include work on the detailed function of these candidate genes and regulatory control of the expression of those genes upon explosion to tumor viruses like MDV.

Materials and Methods
Experimental design. One-day old chickens from two highly inbred lines were sampled for an MDV challenge trial in this study. One genetic line is known as line 6 3 and the other, as line 7 2 . The two genetic lines were developed and have been maintained at the USDA, Agricultural Research Service, Avian Disease and Oncology Laboratory (ADOL) in East Lansing, Michigan. The lines 6 3 and 7 2 share a common major histocompatibility complex (B*2) haplotype but are resistant and highly susceptible to MD, respectively 15 . On the day of hatch, chickens from each line were randomly divided into MDV challenge group and control group. Each of the chicks in the MDV challenge groups of both lines was inoculated intraabdominally with 500 plaque-forming units of 648A passage 10 MDV at day 5 post hatch. No inoculation was implemented in the control groups. Three chickens from each group were randomly euthanized at 10 (latency period) and 21 dpi (reactivation period), respectively. Spleen samples were individually collected, immediately placed into RNAlater solution (Qiagen, Valencia, CA, USA), and stored at −20 °C until total RNA extraction.
All chickens used in this study were housed in a BSL-2 experimental facility during the trial. Feed and water were supplied ad libitum. The chickens were observed daily throughout the entire duration of the experiment. The animal challenge experiment was approved by the USDA, Avian Disease and Oncology Laboratory Institutional Animal Care and Use Committee (IACUC). The IACUC guidelines established and approved by the ADOL IACUC (April 2005) and the Guide for the care and use of Laboratory Animals by Institute for Laboratory Animal Research (2011) were closely followed throughout the experiment. RNA sequencing. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. RNA concentration and quality were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Equal amount of RNA samples from three biological replicates within each line each treatment group were pooled in preparation to construct standard cDNA libraries using Illumina TruSeq kits and reagents following the manufacturer's protocol for deep sequencing. The libraries were sequenced on an Illumina HiSeq2000 sequencer for single end 50 base sequencing run. The post sequencing processes, including image analysis, base calling, and Q-Score calculation, were carried out using Real Time Analysis (v1.13.48); read demultiplexing and conversion to final FASTQ files, using CASAVA (v1.8.2) software tools (Illumina Inc., San Diego, CA, USA). The library preparation, RNA sequence read extraction, and preliminary read quality control were performed at the Research Technology Support Facility, Michigan State University.
Mapping and gene expression quantitation. Sequence adaptors were removed in the first quality control process using Trimmomatic (version 0.32) software 65 to obtain the pass-filter (PF) reads. Low quality bases were further trimmed from the PF reads using custom Python scripts eliminating the first 15 nucleotides. Sickle (v1.33) 66 was used with a sliding window average score of 30 in removing reads with "N"s, and minimum read length of 30 bps, and ended with clean reads. The clean reads were then used to map to the chicken genome reference (galGal4) using TopHat2 (v2.0.12) 67 and Bowtie2 (v2.2.3) 68 with default parameters. Transcript abundance and differential expression of genes were estimated with Cufflink (v2.2.1) 69 . FPKM values were obtained to quantify relative expression of transcripts.
Analyses of DEGs between treatment and chicken line groups. The number of reads per gene for each sample were counted using HTSeq 70 . In each of the pairwise comparisons (between infected birds and control in each line, and between the two lines with and without MDV infection), DEGs were identified by using the DESeq R package (2.1.0) and selected using a filter criterion of FDR <0.05 and FC >2. To better understand the functional involvements of these DEGs, g:Proflier (http://biit.cs.ut.ee/gprofiler/index.cgi) 71 were used for the gene annotation, GO term and pathway enrichment analyses.
Droplet Digital tM PCR validation of gene expression. To validate sequencing data, three genes from each of the treatment groups were selected and re-evaluated on a Droplet DigitalTM PCR (QX200TM ddPCR system; Bio-Rad Laboratories, Inc., Hercules, CA, USA). The primers for ddPCR were designed with Primer3Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi/) and are listed in Table S10. A total of seven genes was selected. The cDNA samples used in ddPCR validation were reversely transcribed from individual RNA samples using the iScript TM RT Supermix Kit (Cat No. 170-8841) and following the manufacturer's instructions (Rio-Rad). Those total RNA samples were the same samples used in cDNA library preparation for the RNA-Seq