Novel Mechanisms Revealed in the Trachea Transcriptome of Resistant and Susceptible Chicken Lines following Infection with Newcastle Disease Virus

ABSTRACT Newcastle disease virus (NDV) has a devastating impact on poultry production in developing countries. This study examined the transcriptome of tracheal epithelial cells from two inbred chicken lines that differ in NDV susceptibility after challenge with a high-titer inoculum of lentogenic NDV. The Fayoumi line had a significantly lower NDV load postchallenge than the Leghorn line, demonstrating the Fayoumi line's classification as a relatively NDV-resistant breed. Examination of the trachea transcriptome showed a large increase in immune cell infiltration in the trachea in both lines at all times postinfection. The pathways conserved across lines and at all three time points postinfection included iCOS-iCOSL signaling in T helper cells, NF-κB signaling, the role of nuclear factor of activated T cells in the regulation of the immune response, calcium-induced T lymphocyte apoptosis, phospholipase C signaling, and CD28 signaling in T helper cells. Although shared pathways were seen in the Fayoumi and Leghorn lines, each line showed unique responses as well. The downregulation of collagen and the activation of eukaryotic translation initiation factor 2 signaling in the Fayoumis relative to the Leghorns at 2 days postinfection may contribute to the resistance phenotype seen in the Fayoumis. This study provides a further understanding of host-pathogen interactions which could improve vaccine efficacy and, in combination with genome-wide association studies, has the potential to advance strategies for breeding chickens with enhanced resistance to NDV.

during outbreaks. Hitchner (4) refers to NDV as a "sleeping giant" that requires close monitoring because the large range of clinical signs of infection with NDV and the genomic diversity of NDV make NDV detection and diagnosis challenging (5,6), the worldwide distribution increases the potential danger of NDV (5), and, although they are uncommon, viruses of low virulence have the ability to mutate and increase in virulence (7).
The range of clinical signs in chickens infected with NDV varies depending on the route of infection, environmental factors, host immunity, and the strain of the virus (8). The strains of NDV are divided into three major categories ranging from low to high virulence: lentogenic, mesogenic, and velogenic (8). The intracerebral pathogenicity index (ICPI) is used to differentiate and classify NDV according to its virulence (8), which is largely determined by the F protein cleavage site (9). Although the strains differ in their virulence, all belong to the same serotype, thus enabling the use of lentogenic viruses as vaccines for protection against the more virulent strains. Although infections by lentogenic strains are rare, clinical signs associated with infections by lentogenic strains are nonspecific and include ruffled feathers, anorexia, decreased egg production, and respiratory infection in birds with a compromised immune system (5,6,10).
Understanding the host response to NDV is necessary in order to generate improved solutions to combat this devastating disease. Examination of host mRNA expression after challenge provides insight into the host-pathogen interaction. Previous studies have compared gene expression levels in response to infection with multiple NDV strains and in host tissues measured from quantitative PCR (qPCR) and microarray data and reported different cytokine expression levels dependent on the strain causing the infection (11)(12)(13)(14)(15). Gene expression analyses performed on tracheal epithelial cells following NDV infection in vitro showed a strong innate immune response (16). The trachea is one of the first tissues that NDV encounters when it is transmitted by the aerosol route and a critical site of the host-pathogen interaction. All strains of NDV are able to replicate in the epithelial cells of the trachea (17), and challenges with either lentogenic or velogenic strains resulted in comparable viral titers in the trachea (18).
Although lentogenic NDV generally does not cause severe clinical signs, gene expression changes due to inflammation, tissue destruction, cell proliferation, and tissue remodeling are expected in the trachea after intranasal inoculation (19). Aerosol delivery can result in deciliation, congestion, goblet cell hyperplasia, edema, and infiltration of heterophils, lymphocytes, and plasma cells in the tracheal mucosa (20,21). If one is seeking mechanisms of host resistance, the trachea is an ideal tissue to examine, as this is one of the initial sites of interaction with the virus.
Host genetics may play a large role in the host-pathogen interaction. Previous studies have shown that resistance to NDV has a genetic basis (22)(23)(24)(25). In the current study, two diverse, inbred chicken lines, the Fayoumi and Leghorn lines, were utilized to identify genetic mechanisms of resistance to NDV. Several studies have compared the immune responses of these two lines. Chickens of the Fayoumi line were found to be relatively resistant to Salmonella, Eimeria, Marek's disease, and avian influenza virus (AIV) compared to those of the Leghorn line (26)(27)(28)(29). Bone marrow-derived dendritic cells from these two lines showed increased phagocytic ability, nitric oxide production, and major histocompatibility complex (MHC) class II surface expression (30). Here the Fayoumi and Leghorn lines were used as a discovery platform to identify mechanisms of resistance to NDV. We hypothesized that the Fayoumis are more resistant to NDV, as shown by a lower viral load and a higher antibody titer. The term "resistance," as used in the current study, is defined as the ability of the host to interfere with the pathogen life cycle (31). The term resistance in this context is not absolute; therefore, the host may still become infected and succumb to infections caused by virulent NDV strains. Comparison of the two lines' responses to NDV may designate genes/pathways associated with the resistance of the Fayoumi line. The information generated by this study will be beneficial for vaccine development and other control strategies.

RESULTS
The Fayoumi line's resistance phenotype was upheld by the NDV viral load, antibody quantification measurements, and sequence data. To examine the effects of an NDV challenge on the two inbred lines with different resistance, the challenged Fayoumi and Leghorn chickens were inoculated with the La Sota NDV strain, while the nonchallenged chickens were given saline solution as a mock infection. Lachrymal fluid was collected from each chicken prior to infection, at 2 days postinfection (dpi), and at 6 dpi for viral quantification. No NDV was detected in any of the birds prior to challenge or in any of the nonchallenged birds at all times (data not shown). In the challenged birds, the viral load significantly decreased from 2 to 6 dpi (P Ͻ 0.0001), and the line had a significant effect on the viral load (P ϭ 0.045) (Fig. 1). At 6 dpi, the Fayoumis had significantly less virus than the Leghorns (P ϭ 0.0122), suggesting that the Fayoumis may have cleared the virus more quickly (Fig. 1). There was no correlation (r ϭ Ϫ0.0097) between each individual's viral load at 2 and 6 dpi, which was in agreement with the findings of previous studies (19).
Serum from blood collected at 10 dpi was used to quantify NDV-specific antibody via enzyme-linked immunosorbent assay (ELISA). The ELISA sample-to-positive (S/P) ratios at 10 dpi were significantly different (P ϭ 0.0007) between the challenged and nonchallenged birds. Due to the large variation in NDV antibody levels, there was no The viral load is shown as the least-square means of the log copy number, measured by qPCR, of the Fayoumis (white) and Leghorns (gray). Error bars represent standard deviations. A Student t test was used to derive the connecting letters report, in which the values for bars that are not labeled with the same letter are significantly different (P Ͻ 0.05). At 2 dpi, data are for 21 Fayoumi and 28 Leghorn chickens, and at 6 dpi, data are for 12 Fayoumi and 20 Leghorn chickens (total n ϭ 81). (B) The antibody titer is displayed as the least-square means of the S/P ratio in the Fayoumis (white) and Leghorns (gray) at 10 dpi, and the error bars represent standard deviations. The connecting letters report was generated by Student's t test. At 10 dpi, data are for 8 challenged Fayoumi chickens, 13 challenged Leghorn chickens, 6 nonchallenged Fayoumi chickens, and 8 nonchallenged Leghorn chickens (total n ϭ 35). significant difference between the lines within the same treatment category. However, the challenged Fayoumis produced more antibodies, on average (P ϭ 0.367), than the challenged Leghorns (Fig. 1).
Sequence reads that did not map to the chicken genome were analyzed to determine if any mapped to genes in the NDV genome (Fig. 2). Viral transcripts were detectable in the challenged birds only at 2 dpi. The main effect of line had a statistically significant impact (P ϭ 0.0264) on the viral transcripts of the challenged Fayoumis and Leghorns at 2 dpi (Fig. 2). As expected, the counts per million (cpm) appeared to be higher for genes at the beginning of the virus genome in both lines.
Effects of NDV challenge on the trachea transcriptome. Epithelial cells were peeled from the trachea of challenged and nonchallenged birds from both lines at all three times and used for RNA isolation, cDNA library construction, and sequencing. The tracheal epithelium sequencing results were similar across treatment groups, and therefore, we did not expect biases due to differences in the number of reads, mapping percentage, or transcriptome coverage percentage among the treatment groups (Table  1). On average, approximately 2 million reads were removed from each sample in the filtering process, nearly 90% of the filtered reads were mappable, and about 75% of the Galgal4 transcripts had at least one mapped read (Table 1).
To compare the host responses to NDV, the responses of challenged and nonchallenged birds within each line were compared at each time point in the experiment.  These comparisons are referred to as the six major contrasts. As time progressed, there were fewer differentially expressed (DE) genes (DEG) between the challenged and nonchallenged birds in both lines (Fig. 3). Overall, more genes were upregulated than downregulated (Fig. 3). Differentially expressed genes that were shared between lines are likely crucial to the response to NDV for all chickens, whereas genes uniquely identified as differentially expressed in the Leghorn or the Fayoumi chickens may be related to the susceptibility and resistance, respectively, of the chickens. The gene for protein kinase C beta (PRKC␤) was the only gene that was uniquely identified to be DE only in the Fayoumis at all three time points. The large difference in the numbers of DEG at 6 dpi in the Fayoumi and Leghorn lines (Fig. 3) suggests that the Fayoumi chickens recovered more quickly from the virus infection, resulting in fewer differences between the challenged and nonchallenged birds.
Challenge with NDV stimulated the infiltration and migration of immune cells. The DEG from the challenged versus nonchallenged birds at each time point and within each line were used for cell type enrichment analysis. At 2 dpi, the top cell types enriched were similar for the two lines (Fig. 4); whole-blood cells were the most enriched cell type for both lines. At 6 dpi, the lines were the most dissimilar, which corresponds to the difference in DEG numbers ( Fig. 3 and 4). The enrichment of cancer-type cell lines at 6 dpi may have been a result of high levels of cell proliferation. Overall, the results suggested a conserved response between the two lines and a strong enrichment of immune-related cells after challenge at all time points but no clear shifts from innate to adaptive immune cell types over time.
A major pathway predicted by Ingenuity pathway analysis (IPA) to be activated in all six major contrasts was the leukocyte extravasation signaling pathway. Figure 5A shows the genes in this pathway that were significantly up-or downregulated within each contrast. On the basis of the expression levels of the genes in this pathway (Fig. 5B), the molecule activity prediction function of IPA predicted, for each contrast, activation of cell mobility and activation of actin cytoskeleton contraction; the IPA output placed more confidence in the prediction of the latter in the Fayoumis than the Leghorns at all time points (not shown).
Temporal changes influenced the pathways that were activated or inhibited in each line. The top canonical pathways predicted by IPA are shown in Fig. 6 for each of the six major contrasts. The six major contrasts were clustered by similarity, which was determined by the resemblance of the activation or inhibition scores of all the pathways. The clustering grouped the major contrasts by dpi, with the Fayoumis at 6 dpi being the outlier ( Fig. 6). At 6 dpi there were only 119 DEG for the Fayoumis to populate the pathways, making it more difficult to obtain significance for activation or inhibition of a pathway. Several key pathways were significantly activated at all time points in both lines: iCOS-iCOSL signaling in T helper cells, NF-B signaling, the role of the nuclear factor of activated T cells in the regulation of the immune response, calciuminduced T lymphocyte apoptosis, phospholipase C signaling, and CD28 signaling in T helper cells. These conserved pathways may be inherently crucial to the chicken's response to NDV in the trachea. The stress-activated protein kinase (SAPK)/Jun N-terminal protein kinase (JNK) signaling pathway was activated only in the Leghorns and was activated at all three time points. This pathway may be associated with the Leghorns' susceptibility to NDV (Fig. 6).
Collagen downregulation at 2 dpi was identified only in the Fayoumi chickens. The network shown in Fig. 7 shows the overwhelming downregulation of collagen, which was predicted by IPA to be associated with various diseases and functions, such as connective tissue disorders, dermatological diseases and conditions, and organismal injury and abnormalities. This network had 25 focus modules and had the secondhighest score of the networks generated from this contrast. Of the 30 genes in this network, only 8 were also DE in the Leghorn chickens (Fig. 7). It does not appear that the downregulation of collagen was limited to a specific collagen type.
Direct comparison of the Fayoumi and Leghorn chickens at 2 dpi gave more insight into the pathways and genes associated with resistance. The Fayoumis and Leghorns had a high degree of similarity in the cell populations enriched at 2 dpi ( Fig.  4), allowing for direct comparison between the two lines. The number of DEG between the Fayoumis and Leghorns was higher in the challenged birds than the nonchallenged birds ( Fig. 8). Myosin-and troponin-related genes, as well as MHC class I genes, were the most upregulated genes in the Fayoumis relative to the Leghorns at this time. The eukaryotic translation initiation factor 2 (eIF2) signaling pathway had the highest positive Z-score (Fig. 8), suggesting strong activation of this pathway in the Fayoumi chickens relative to the Leghorn chickens.
Examining the genes that were impacted by challenge and line offered potential candidate genes of interest. A contrast was written to compare the interaction between challenge and line at each time point. At each time point, only DEG from the challenged versus nonchallenged contrast within each line were examined for the interaction (line · challenge). The majority of the genes with a false discovery rate (FDR) of Ͻ0.05 for the interaction were from 6 dpi, and nearly half of all genes were uncharacterized proteins (Table 2). Further investigation into the function of these uncharacterized proteins would be useful to provide more insight into the resistant and susceptible phenotypes of Fayoumi and Leghorn chickens.
The RNA-seq results were validated with an independent test. The Fluidigm Biomark HD system was used as a method of high-throughput quantitative PCR to serve as validation for the transcriptome sequencing (RNA-seq) results. The expression of 44 genes (Table 3) was analyzed. Figure 9 shows a correlation of 0.91 between the FIG 5 Leukocyte extravasation signaling pathway. For all genes, the absolute LFC was Ͼ1 and the FDR was Ͻ0.05. (A) Gene expression heat map of genes in the leukocyte extravasation signaling pathway generated by contrasting the challenged and nonchallenged birds within each time and line. The contrasts, listed atop each column, were clustered by similarity on the basis of genes (blue text) that were upregulated (red), that were downregulated (green), and that had no change (gray). Genes that were not DE in any of the comparisons were removed from the graph. This pathway was activated in all contrasts (orange), as predicted by IPA. (B) The leukocyte extravasation signaling pathway in the Fayoumis at 2 dpi. Genes were upregulated (red), downregulated (green), or not DE (gray) in the challenged relative to the nonchallenged chickens. Genes in white were not found in the data set because they were removed during normalization, unrecognized by IPA, or not found in the chicken genome. The shapes of the genes represent their function. The genes within a parent node that were DE in any of the six contrasts are shown individually next to the parental node and are denoted with smaller molecules and labels.
two methods across all six major contrasts. The high correlation serves as a validation of the RNA-seq technology applied in this study.

DISCUSSION
Examining changes in the transcriptome of the tracheal epithelia between challenged and nonchallenged birds offered valuable insight into the host response to NDV. This study predicted the enrichment of immune cells in the trachea in challenged birds of both lines at all times. One pathway that was instrumental in leukocyte migration was leukocyte extravasation signaling. On the basis of the expression of genes within this pathway, IPA predicted the activation of actin cytoskeleton contraction and cell mobility. A previous study examining gene expression levels in the lung after challenge with AIV in the same experimental lines showed an upregulation of FIG 6 Pathways predicted to be activated or inhibited in response to NDV in the six major contrasts. For all genes, the absolute LFC was Ͼ1 and the FDR was Ͻ0.05. (A) The heat map, generated by Qiagen IPA software, shows pathways that were either activated (orange), inhibited (blue), or unchanged (white) when the findings for the challenged and nonchallenged birds within each line were contrasted at the three separate time points. The pathways are sorted by hierarchical clustering, and the contrasts (columns) are clustered by similarity. Pathways are filtered by a Ϫlog(P value) of Ͼ2.5 and a Z-score of Ͼ1. The complete names of the pathways for entries with ellipses are provided in panel B. (B) The top pathways [those with a Ϫlog(P value) of Ͼ2.5 and a Z-score of Ͼ2] for each contrast are shown. The total number of genes in each pathway is found at the top of each bar, which is divided into genes that are downregulated, upregulated, and not DE in the gene data set according to the percentage of the total number of genes in each pathway. The orange lines represent the Ϫlog(P value) of each pathway. NFAT, nuclear factor of activated T cells; SAPK, stress-activated protein kinase; JNK, Jun N-terminal protein kinase; PI3K, phosphatidylinositol 3-kinase; iNOS, inducible nitric oxide synthase; MAPK, mitogen-activated protein kinase; PDGF, platelet-derived growth factor; IL-9, interleukin-9; GM-CSF, granulocyte-macrophage colony-stimulating factor. actin filament-based movement in the Fayoumis (29). Actin cytoskeleton activation is crucial for cell migration, as leukocytes use actin polymerization and actomyosin contraction for movement. Notably, paramyxoviruses like NDV also require actin to replicate efficiently (32,33). If the migrating cells become infected with the virus, a conundrum may be created for the host. In one study, there was a downregulation of actin-related genes after challenge with NDV in vitro (34). Due to challenge with infectious bronchitis virus (IBV), the proteins involved in cytoskeleton organization differed in expression in the trachea depending on IBV virulence (35). The host must find a way to balance the costs and benefits of cytoskeleton upregulation with defense against viral infection. The regulation of cytoskeleton-and actin-related genes may be time dependent and may be a critical target for future study.
Host-pathogen interactions at the site of infection can determine the severity of both the pathology and the overall infection. The use of two unique inbred lines that differ in their relative resistance offered a valuable tool for comparison. Fayoumis had significantly less detectable virus at 2 dpi (measured by RNA-seq) and 6 dpi (measured by qPCR), which may result in faster viral clearing and which, overall, supports their relatively resistant phenotype compared to the phenotype of the Leghorns. A decreased viral load corresponds to less viral shedding and, therefore, fewer viral particles in the environment and a lower transmission potential (36). In a separate study with the same two genetic lines, although all birds succumbed to infection with velogenic NDV, the mean survival time was longer for the Fayoumis than the Leghorns at all four doses tested and was significantly longer at the two higher doses (P. Miller and C. Afonso, SEPRL, Athens, GA, personal communication, 2016). As the time postinfection increased, the numbers of DEG between challenged and nonchallenged birds decreased, suggesting that the chickens were recovering. The Fayoumis had many fewer DEG at 6 dpi than the Leghorns, which may be further evidence that the Fayoumis recovered more quickly. Although most pathways had similar predicted levels of activation or inhibition in the two lines, the lines often differed in gene expression levels within the pathway (Fig. 5). The current study identified several genes whose expression levels were impacted significantly by both line and challenge (Table 2), and some of these overlapped those from a previous   (29). These genes may be of particular importance in the host response to both AIV and NDV. The wide range of pathogens (parasites, bacteria, viruses) to which the Fayoumis have shown relative resistance suggests that the mechanism(s) for their resistance may function early in the immune response (26,27,29). Downregulation of collagen at 2 dpi in the challenged Fayoumis could be partially responsible for their resistant phenotype. One study found lower levels of expression of collagen mRNA in the trachea of resistant chickens (noninbred Cornell chickens) than susceptible chickens (inbred White Leghorn chickens) after challenge with IBV (37). The downregulation of collagen may impact apoptosis (38), immune cell migration (39)(40)(41), and T-cell activation (42)(43)(44)(45)(46); impacting these processes would likely impact the virus and its pathogenesis.
A direct comparison of the challenged Fayoumis and Leghorns at 2 dpi also gave insightful results. Integrin signaling was more highly activated in the Leghorns than the a The primers were used/analyzed previously (62).
Fayoumis (Fig. 8). This may be directly related to the downregulation of collagen and its relationship with integrin (38). Laminin, another key extracellular matrix protein (39), was also upregulated in the Leghorns relative to the Fayoumis. The impact of collagen was shown in two separate contrasts. The downregulation of collagen between challenged and nonchallenged Fayoumis and the activation of collagen-related pathways in the Leghorns relative to the Fayoumis give more confidence in the importance of collagen at 2 dpi. The activation of the protein kinase R (PKR)/eIF2 signaling cascade is known to inhibit NDV replication (47), likely due to the actions of the antiviral protein PKR. This pathway, involved in mRNA translation, was more upregulated in the Fayoumis than the Leghorns (Fig. 8). This finding was in agreement with the Fayoumis' resistant phenotype and may be related to the lower viral load seen at 6 dpi. The phosphorylation of eIF2␣ results in apoptosis due to a global termination of transcription (48). Apoptosis-related pathways were activated, including T lymphocyte apoptosis, p53 signaling, retinoic acid-mediated apoptosis signaling, and apoptosis signaling (Fig. 6). The eIF2 signaling pathway was also identified to be differentially methylated between chickens that were either vaccinated or not for infectious laryngotracheitis (48). It is clear that this pathway is critical in the host defense against multiple pathogens.
In a previous study, chicken stem cells underwent rounds of selection with high titers of the La Sota strain; the levels of expression of these selected cells led to increased levels of apoptosis and decreased neuroactive ligand-receptor interaction compared to those in the unselected cells (49). The Fayoumis showed higher levels of activation of apoptosis signaling at 2 dpi (Fig. 6), and when the lines were directly compared, the Leghorns showed more activation of the neuregulin signaling pathway at 2 dpi (Fig. 8). The overlap in these pathways suggests that increased levels of apoptosis and decreased neurologically related pathways may be important mechanisms of host resistance to La Sota.
The trachea transcriptome exposed potential key mechanisms for the host defense against NDV. As expected, there appeared to be a large increase in immune cell infiltration in the trachea in both lines at all times examined. Although shared mechanisms were seen in both the Fayoumis and Leghorns, the lines showed unique responses as well. This study has demonstrated the Fayoumis' classification as being FIG 9 RNA-seq results validated using the Fluidigm Biomark HD system. This scatter plot shows the agreement between the LFC measured by RNA-seq and the LFC measured by the Fluidigm Biomark HD system. The LFCs were generated from the six major contrasts for the challenged and the nonchallenged birds at each time and in each line for 43 genes in total. Each dot represents the LFC of a gene in a particular contrast; there are 258 points in total. The pairwise correlation between the two methods was 0.91. H6PD was the housekeeping gene used for normalization of the data obtained with the Fluidigm Biomark HD system. See Table 3 for primer information.
relatively resistant to NDV compared to the Leghorns, as determined by viral load and viral sequence counts and suggested by antibody titers and DEG numbers. This study identified possible mechanisms for the Fayoumis' resistance, including the downregulation of collagen and the increased activation of eIF2 signaling, which can be used to improve vaccine development and highlight genes for which beneficial genetic variance can be used to inform breeding decisions. Cytokine expression levels are time, tissue, and strain dependent (11)(12)(13). It is critical that future breeding decisions be made from data from studies with multiple tissues, host genetic backgrounds, NDV strains, and time points. The current study is a first step in understanding and identifying possible mechanisms of resistance to NDV. Future studies to examine the gene expression of Fayoumis and Leghorns challenged with velogenic NDV are necessary to better determine the applicability of the results presented here.

MATERIALS AND METHODS
Ethics statement and animals used. The methods used in the present study were approved by the Iowa State University IACUC (log number 1-13-7490-G). Fayoumi (M 15.2) and Leghorn (GHs 6) chicken lines from Iowa State University Poultry Farm (Ames, IA) were used in this study. The Fayoumis and Leghorns each have an inbreeding coefficient of Ͼ99.9%, as determined by microsatellites (50), and 99.95%, as determined by resequencing (51). The Fayoumi line originated from the Fayoum region of Egypt, whereas the Leghorn line is representative of chickens from the U.S. egg-laying industry of the 1950s. The chickens were given ad libitum access to water and feed throughout the study. Supplemental heat at 35°C was provided to the birds starting on the day of hatch, and the temperature was stepped down 2 to 3°C every few days until it reached 24°C when the chicks were 29 days of age and remained at 24°C for the remainder of the study. The light cycle began at 23 h of light and 1 h of dark on the day of hatch, and the number of hours of light was gradually reduced until it reached 13.5 h of light when the chicks were 24 days of age and remained at that amount for the remainder of the study.
Virus. A live attenuated Newcastle disease vaccine, the type B1 La Sota lentogenic strain, was propagated once in 10-day-old embryonated specific-pathogen-free (SPF) chicken eggs inoculated via the allantoic sac (52). The allantoic fluid was harvested after 2 days, tested for virus presence by determination of hemagglutinating activity (HA), and centrifuged at 8,000 rpm for 30 min to clear the debris. Virus titration was performed in SPF embryonated eggs as previously described (53). The final titer of the undiluted virus suspension was 10 8 50% embryo infectious dose (EID 50 ) per ml. The virus suspension was aliquoted and stored at Ϫ80°C until its use.
Experimental design. At the time of hatch, straight run chickens of both lines were randomly assigned to one of four rooms in a biosafety level II facility. At 21 days of age, half of the chickens (two rooms, n ϭ 49) were inoculated with 200 l of 10 7 EID 50 of La Sota NDV, with 50 l being inoculated into each of the eyes and nostrils. The nonchallenged group (2 rooms, n ϭ 40) was given 200 l phosphatebuffered saline (PBS), via the same route, as a mock infection. Chickens that received the La Sota virus are referred to throughout this article as "challenged," and those that received the mock infection are referred to as "nonchallenged." To collect lachrymal fluid prechallenge (n ϭ 89) and at times of 2 (n ϭ 89) and 6 (n ϭ 62) dpi, fine, crystalline sodium chloride was placed on the eye of each chick to collect lachrymal fluid with a pipette. To quantify antibody levels, venous blood samples were collected prechallenge and at 10 dpi. At 2, 6, and 10 dpi, 24, 30, and 35 chicks, respectively, were euthanized with sodium pentobarbital solution for tissue harvest. The entire trachea was harvested from all of the challenged and nonchallenged birds sacrificed at each dpi and placed into RNAlater solution (Thermo Fisher Scientific, Waltham, MA) for short-term storage. Within a week, the sheet of epithelial cells from each trachea was carefully removed by peeling it from the interior surface of the trachea with forceps. This epithelial cell lining was placed into a Ϫ80°C freezer for long-term storage, and the remaining tissue was discarded.
NDV antibody. Blood samples collected prechallenge and at 10 dpi were centrifuged at 15 rpm for 5 min to collect the serum supernatant. Serum NDV antibody levels were measured using an Idexx NDV ELISA for chickens (Idexx Laboratories, Inc., Westbrook, ME). The sample-to-positive (S/P) ratio was calculated by subtracting the negative-control value from the average absorption value for each sample and dividing that value by the positive-control value minus the negative-control value. The positive and negative controls were supplied with the Idexx kit. The S/P ratio was calculated from the average absorbance value for each sample, which was run in duplicate. A standard least-squares, effect leverage analysis in JMP (JMP Group Inc., San Francisco, CA) including the main effects of line, treatment, and line · treatment was used; a Student t test generated the connecting letters report. Both tests are parametric. Viral load. Viral RNA was isolated from the chicken lachrymal fluid collected prechallenge and at 2 and 6 dpi using a MagMAX-96 viral RNA isolation kit (Life Technologies, Carlsbad, CA) and was quantified using an LSI VetMAX Newcastle disease virus real-time PCR kit (Life Technologies, Carlsbad, CA) on an MJ Research Opticon2 qPCR machine (Bio-Rad, Hercules, CA). The primers were directed against the matrix protein (M) gene, and dilutions of inoculum virus were used as standards. The log copy number per 1 l of isolated viral RNA was calculated for the mean for each sample, which was run in duplicate. Using JMP statistical software (JMP Group Inc., San Francisco, CA), a standard least-squares, effect leverage test including the main effects of line, dpi, and line · dpi was performed to determine significance, and a Student t test was utilized to generate a connecting letters report. The log copy number of the NDV genome does not show the number of infectious particles.
RNA isolation and cDNA library construction. RNA was isolated from the tracheal epithelial cells using an RNAqueous kit (Thermo Fisher Scientific, Waltham, MA), and samples were then treated with DNase using a DNA-free kit (Thermo Fisher Scientific, Waltham, MA). The quality of all samples was validated (RNA quality number Ͼ 8.0) using the Fragment Analyzer (Advanced Analytical Technologies, Ankeny, IA). A 500-ng input for each sample was utilized to construct the cDNA library using the high-throughput protocol in the TruSeq RNA sample preparation guide (v2; Illumina, San Diego, CA). Following cDNA library validation using the Fragment Analyzer (Advanced Analytical Technologies, Ankeny, IA), samples were sequenced on a HiSeq2500 platform (Illumina, San Diego, CA) for 100-bp, single-end reads (DNA Facility, Iowa State University, Ames, IA).
RNA-seq design. This study included 3 main factors: treatment (challenged, nonchallenged), line (Leghorn, Fayoumi), and time (2, 6, 10 dpi), resulting in 12 treatment groups. Each treatment group was represented by 4 biological replicates, except for the Leghorns at 2 dpi, for which there were 3 nonchallenged and 5 challenged chickens included in the analysis. Equal numbers of males and females were used for each treatment group. Samples from birds in the 12 treatment groups were balanced across four lanes of the flow cell and randomly assigned an index.
Data analyses. The Discovery Environment of iPlant Collaborative (54) was utilized for data processing. First, the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was utilized for an initial quality check and assessment of all sequence data. The Illumina TruSeq adapter was recognized as an overrepresented sequence in some samples. For consistency, the adapter sequence and the individual multiplexing index were removed from all samples using the FASTX Clipper program. In addition, FASTX removed any sequence less than 30 bp in length and filtered each read such that 80% of all base pairs within a read had a Phred quality score of 30 or higher. The remaining high-quality reads were used as input for the TopHat2 aligner (55) for alignment to the Galgal4 reference genome (assembly GCA_000002315.2). The Galgal4 (assembly GCA_000002315.2) and GTF files were downloaded from the Ensembl genome browser. In TopHat2 (55), default settings were used, except the minimum isoform fraction was lowered to 0.10 (the default is 0.15) to identify more alternatively spliced transcripts. The HTSeq library (56) was used to count the reads aligned to each transcript (the intersection-nonempty mode was used). Statistical analysis of the count data was performed using the generalized linear model option in the edgeR package (57). The main effects of line, treatment, and time were included in the model. The contrasts comparing the challenged and nonchallenged birds within each line at each time are referred to as the six major contrasts (Fayoumi, 2 dpi, challenged versus nonchallenged; Leghorn, 2 dpi, challenged versus nonchallenged; Fayoumi, 6 dpi, challenged versus nonchallenged; Leghorn, 6 dpi, challenged versus nonchallenged; Fayoumi, 10 dpi, challenged versus nonchallenged; Leghorn, 10 dpi, challenged versus nonchallenged).
The differentially expressed genes (DEG; FDR Ͻ 0.05) identified in either chicken line when comparing challenged to nonchallenged included 1,688, 965, and 180 DEG at 2, 6, and 10 dpi, respectively (Fig.  3). For these genes, contrasts were written to determine the significance of the interaction between line and challenge at each time. The P values generated from the interaction contrast were subjected to a Benjamini and Hochberg multiple testing correction (58).
NDV sequence discovery. The NDV La Sota sequence (GenBank accession number JF950510.1) was downloaded from NCBI and used to make a GFF file. Using the Burrows-Wheeler alignment (BWA) (59) under default settings, the unmapped reads were aligned to the NDV genome (GenBank accession number JF950510.1). The HTSeq library (56) was used to count the number of reads that aligned to each viral gene. The cpm was calculated by taking the number of reads that aligned to each viral gene, dividing by the total number of unmapped reads, and multiplying by 1 million. There was no statistically significant difference between the number of unmapped reads in the challenged Leghorns and Fayoumis at 2 dpi (P ϭ 0.4). The data were imported into JMP statistical software (JMP Group Inc., San Francisco, CA) and analyzed with a least-squares analysis and Student's t test to generate the connecting letters report.
Gene expression analysis. The trachea is a heterogeneous tissue composed of multiple cell populations that likely change under different treatment conditions. To determine cell type enrichment, the transcript identifiers of the DEG (FDR Ͻ 0.05) for each of the six major contrasts were converted into their associated gene names using BioMart software (http://uswest.ensembl.org/biomart/martview/). The six gene lists were then input into the CTen database (http://www.influenza-x.org/~jshoemaker/cten/) (60). The CTen database is based on human and mouse data and has previously been used for chickens (61).
For pathway analysis, Qiagen's Ingenuity Pathway Analysis (IPA; Qiagen, Redwood City, CA) was used. The identifiers of the transcripts unrecognized by IPA were converted to a recognized identifier using BioMart software (http://uswest.ensembl.org/biomart/martview/) to ensure that the maximum number of transcripts was included for pathway analysis. The genes used for pathway analysis had an FDR of Ͻ0.05 and an absolute log 2 fold change (LFC) of Ͼ1, as estimated by the edgeR package. Fluidigm Biomark HD system. To validate the RNA-seq results, the Fluidigm Biomark HD system (Fluidigm, South San Francisco, CA) was used as a method of high-throughput qPCR. The input RNA (25 ng/l) came from the same isolate used for cDNA library construction. The 48 RNA samples were converted into cDNA, preamplified for 12 cycles with all 48 primer pairs (Table 3), treated with exonuclease, diluted 10-fold, and analyzed on the 48.48 integrated fluidic circuit (IFC) for a total of 2,304 sample · gene tests. The genes were chosen to represent the range of log fold changes; not all genes were DE. Ribosomal protein L4 (RPL4) and hexose-6-phosphate dehydrogenase (H6PD) were two housekeeping genes with a high raw cycle threshold (C T ) value correlation (r ϭ 0.9). Since its C T values were closer to the average C T values for all primers and samples, H6PD was used as the housekeeping gene for the analysis. Data were filtered by Fluidigm quality assessment and melting curve consistency. If, after these filtering steps, primer pairs did not have more than 28 samples remaining, all measurements for that gene were excluded. A total of 9.5% of the data, including those for four genes (MZB1, 5_8S_rRNA, UbI, and GAL2), were excluded from the analysis. Table 3 shows the genes and primers analyzed in this study.
For all six major contrasts, the LFC in expression for each gene was calculated by the Ϫ2 ϪΔΔCT method. The correlation between the RNA-seq-generated LFC and the Fluidigm Biomark HD systemgenerated LFC was determined with a pairwise correlation across all genes and all six major contrasts for 258 pairwise comparisons in total.