Novel analysis of the Harderian gland transcriptome response to Newcastle disease virus in two inbred chicken lines

Behind each eye of the chicken resides a unique lymph tissue, the Harderian gland, for which RNA sequencing (RNA-seq) analysis is novel. We characterized the response of this tissue to Newcastle disease virus (NDV) in two inbred lines with different susceptibility to NDV across three time points. Three-week-old relatively resistant (Fayoumi) and relatively susceptible (Leghorn) birds were inoculated with a high-titered (107EID50) La Sota strain of NDV via an oculonasal route. At 2, 6, and 10 days post infection (dpi) Harderian glands were collected and analyzed via RNA-seq. The Fayoumi had significantly more detectable viral transcripts in the Harderian gland at 2 dpi than the Leghorn, but cleared the virus by 6 dpi. At all three time points, few genes were declared differentially expressed (DE) between the challenged and nonchallenged birds, except for the Leghorns at 6 dpi, and these DE genes were predicted to activate an adaptive immune response. Relative to the Leghorn, the Fayoumi was predicted to activate more immune pathways in both challenged and nonchallenged birds suggesting a more elevated immune system in the Fayoumis under homeostatic conditions. Overall, this study helped characterize the function of this important tissue and its response to NDV.

The Harderian gland is a small tissue located behind the eyes of the chicken. In a previous study of adult chickens, the average weight was found to be 84.4 mg 1 . The major functions of this tissue have not been fully characterized but are thought to involve lubrication of the nictitating membrane 2 and the local immune response, serving as a site of antibody production 3 . The cells of the Harderian gland are both bursa-and thymus-dependent 4,5 . Large numbers of heterophils and plasma cells accumulate in the chicken Harderian gland by two weeks of age 2 , and this is a location of terminal B cell maturation 6 . Antigenic stimulation of the Harderian gland also increases the number of plasma cells 7,8 . Although some studies have shown evidence of low numbers of T cells residing in the Harderian gland 9,10 , others saw an increase in T cell numbers after infection with Newcastle disease virus (NDV) 11,12 . There have been several studies on this unique tissue in the past 50 years, but there has been very little recent characterization of the Harderian gland with new technologies. To the authors' knowledge, the Harderian gland transcriptome has never been reported prior to this study.
NDV is a major global problem in poultry. This virus can cause high levels of mortality when birds are not vaccinated, or vaccinated improperly 13 . Eye-drop vaccination is a common method and it is unknown how the transcriptome of the closest immune tissue responds to vaccination.
Genetics plays a role in the chicken's response to NDV [14][15][16] and to the live virus vaccine [17][18][19] . A previous study has identified two inbred chicken lines that differ in their susceptibility to NDV 17 . The Fayoumi (relatively resistant) had significantly lower viral load in the lachrymal fluid at 6 dpi, numerically higher serum antibodies at 10 dpi, and significantly lower viral transcript counts in the tracheal epithelial cells than the Leghorn (relatively susceptible) after challenge with lentogenic NDV 17 . Studying these lines' response to NDV in the Harderian gland can offer insights into both vaccine development and host resistance. We predict these two lines will respond differently to NDV in the Harderian gland transcriptome in terms of viral transcript counts, numbers The principal component analysis plot shows clear clustering by line. The principal component analysis (PCA) plot generated by pcaExplorer shows that principal component (PC) 2 separated the samples clearly by line and accounted for 8.6% of the variance seen in the Harderian gland transcriptome (Fig. 2). PC1 accounted for a large portion of the variance (51.1%), however, samples did not clearly cluster based on any known parameters. PC1 may be related to dpi because some time points clustered very closely on PC1 (Fig. 2). The differentially expressed genes between the nonchallenged birds over time within each line are included in Supplementary Table S1.
Challenged vs. nonchallenged birds by line and time. The challenge resulted in few DEG at most time points within each line (Table 2). Of the few DEG, several have known immune function including STAT1, Mx1, and TLR3, which were all upregulated due to challenge ( 6 dpi had much larger numbers of DEG (Table 2), which may be related to lack of viral clearance at 6 dpi in the Leghorns as shown by higher viral cpm values (Fig. 1). Pathway analysis was performed to further characterize the response to challenge in the Leghorns at 6 dpi because of the large number of DEG ( Table 2). The top five canonical pathways with a z-score >0.01 generated by IPA included the Th2 Pathway, Th1 Pathway, iCOS-iCOSL Signaling in T Helper Cells, CD28 Signaling in T Helper Cells, and PKCθ Signaling in T lymphocytes. All of these pathways were predicted to be activated and involved the T cell response. The B Cell Development Pathway was also significant and predicted to be activated, but was not among the top 5 canonical pathways. The upstream regulator T cell receptor (TCR) was predicted to be activated due to the increased expression of the genes (Fig. 3). The expression levels of these genes were also predicted to activate quantity of B lymphocytes, formation of lymphoid tissue, and cell proliferation of T lymphocytes and inhibition of replication of RNA virus (Fig. 3). Thus, in response to NDV infection, at 6 dpi, the Leghorns were activating both the cell mediated and humoral arms of their adaptive immune response in this unique tissue.

Genetic line differences were apparent in the Harderian gland transcriptome in both challenged
and nonchallenged birds at all times. The   the challenged Fayoumis and Leghorns was larger than between the age-matched, nonchallenged Fayoumis and Leghorns (Table 3). This suggested that their transcriptomes became more different after challenge with the virus, even though within-line contrasts of challenged vs. nonchallenged birds showed few DEG at this time ( Table 2). The extremely tight clustering of transcriptomes of the challenged Leghorns at 2 dpi ( Fig. 2) may have increased the ability to call genes as differentially expressed between the challenged Fayoumis and Leghorns. Overall, 58 genes were consistently differentially expressed in the Fayoumis and Leghorns at all times and in both challenged and nonchallenged birds. These DEG may be relevant to breed characterization. Somatotropin, for example, was always more highly expressed in the Fayoumis (Table 3) and has previously been associated with viral resistance to Marek's disease in chickens 20 . Further analysis of the DEG between the Fayoumi and Leghorn revealed which pathways were more highly activated within each line (Fig. 4). Many inflammatory pathways were more activated in the challenged Fayoumis at 2 dpi, and this may be related to the amount of viral transcript counts present at that time. Interestingly, most of the pathways that were differentially activated between the two lines appear to be immune related (Fig. 4). Pathways in which the predicted activation was of different direction in challenged and nonchallenged birds are of particular interest.    RNA-seq results were validated by high throughput qPCR. Fluidigm Biomark was used as a method of high throughput qPCR to validate the RNA-seq methods by correlating the log fold change (LFC) calculated by both techniques. Primers were previously published 17 . The six challenged vs. nonchallenged contrasts (within each time and line), were used to compare the RNA-seq LFC with the Biomark LFC calculated via the −2 −ΔΔCT method (Fig. 5). The correlation (r = 0.77) was strong.

Discussion
There was a large number of DEG between the challenged Fayoumis and Leghorns at 2 dpi. This was surprising, because there were few DEG in the contrasts of challenged and nonchallenged birds within-line at 2 dpi. The PCA can explain some of this result. The transcriptome of challenged Leghorns at 2 dpi clustered very tightly together with the two largest principal components (which together account for 60% of the variance), much tighter than the nonchallenged Leghorns at 2 dpi. The challenged Fayoumis clustered more tightly than the nonchallenged Fayoumis as well at 2 dpi. A multi-dimensional-scaling plot showed the same clustering patterns as the PCA plot (data not shown). Less within-group variation allows for the identification of more DEG. The PCA plot also showed the tight clustering of the nonchallenged Leghorns at 6 dpi, which may have allowed for the detection of more DEG when contrasting the challenged and nonchallenged birds at this time. The samples within Leghorn treatment groups tended to diverge with time, whereas in the Fayoumis, the samples clustered more tightly as time progressed across PC1. It is unclear what drives some treatment groups to cluster more tightly than others; PC1 accounts for a large portion of variance and cannot be identified from the study parameters. Previous studies have shown the influx of plasma cells and T cells to the Harderian gland after infection with NDV 5,11,12 . In the Harderian gland, plasma cells have been shown to synthesize mainly IgA 5 , which has been detected as early as 4 dpi in the saliva of Leghorns 21 after infection with NDV. However, no IgA was detected in the lachrymal fluid of the challenged Fayoumis or Leghorns at 2 & 6 dpi (data not shown).
The large number of DEG seen between the challenged and nonchallenged Leghorns at 6 dpi shows the activation of the adaptive immune response in these birds. These DEG predicted the activation of several T cell-related pathways. The activation of T cell-related pathways after infection was also seen in the trachea epithelial cells of these birds 17 , but not in the lung, which had no detectable viral transcript counts 18 . T cell-related pathways appear to be activated at sites of viral replication, which is likely a result of an increased T cell population. In another study, after infection with Hitchner B1 NDV strain, all T cell subsets increased at least 2-fold in the Harderian gland, and CD8+ T cells increased 6-fold 12 .
There might be a relationship between the large number of DEG and the viral transcript counts in the Harderian gland of the Leghorns at 6 dpi. Viral persistence in the Leghorns could cause the large DEG response at 6 dpi in the Leghorns that was absent in the Fayoumis, because the Fayoumis cleared the virus by 6 dpi. It is unknown, however, if the viral transcript counts correlate with the number of infectious viral particles. Surprisingly, the Fayoumis had significantly more detectable viral transcript counts, as estimated from the RNA-seq, at 2 dpi than the susceptible Leghorns in the Harderian gland. This result should be explored further to determine whether greater viral transcript counts is advantageous for mounting an appropriate immune response in the Harderian gland. Overall, the cpm for each viral gene detected in the Harderian gland was lower than that detected in the trachea epithelial cells of these birds 17 . The different cell composition of these tissues may be one of the causes. Likely, a higher percentage of cells were susceptible to lentogenic NDV replication in the trachea epithelial cells than the Harderian gland, as lentogenic NDV requires a protease found in epithelial cells to cleave the Fusion protein 22 . Also, unlike the trachea, in the Harderian gland the viral gene cpm did not decrease from 3′ to 5′ as would be expected 17 . The high cpm of the F and HN genes in the Harderian gland of challenged Fayoumi at 2 dpi may have facilitated the production of more neutralizing antibodies to these external structural proteins.
The two genetic lines clearly responded differently to the virus in the number of DEG and in the amount of detectable viral transcript counts. When the lines were directly compared, several immune related pathways differed in their relative activation. Relative to the Leghorns, the Fayoumis had more immune-related pathways that were in an activated state under resting conditions (nonchallenged), including eNOS signaling, Th1 pathway, LXR/RXR Activation, B cell receptor signaling, FcγRIIB signaling in B lymphocytes, and production of nitric oxide and reactive oxygen species. These pathways represent each facet of the immune system: innate, cell mediated, and humoral immunity. The relative activation of these immune pathways in the nonchallenged Fayoumis may be a reason for the low numbers of DEG found between the challenged and nonchallenged Fayoumis at the time points measured. Differences in immune related pathways of nonchallenged Fayoumis and Leghorns was also observed in the spleen 19 . It is unknown whether the reason for the Fayoumi's resistant or Leghorn's susceptible phenotype is due to their response to pathogen or their constitutive expression levels. The increased relative activation of immune related pathways under homeostatic conditions may be partially responsible for the favorable phenotype observed in the Fayoumis. The Fayoumis may be more equipped to handle a viral infection under resting conditions and, therefore, do not need to mount a large immune response to clear the virus, as assessed by changes in gene transcription.
The genes that were significantly impacted by the challenge*line interaction are of particular interest. Of those genes, some, including RPLP2, DSCAM, and CPA6, have been shown to have immune function. Flaviviruses require RPLP2 for infection, because it is necessary for viral translation 23 . In arthropods, DSCAM plays a crucial The resulting differentially expressed genes with FDR <0.10 were used as input into IPA for pathway analysis. Contrasts from 10 dpi were removed because no pathways were significant. Heat map shows pathways predicted to be relatively more activated in the Fayoumi (orange) and Leghorn (blue). The more intense a square in the heat map, the higher the z-score and prediction confidence. White squares represent contrasts with no or too few DEG within that pathway to make a prediction (z-score = 0).
SCIentIfIC REPORTS | (2018) 8:6558 | DOI:10.1038/s41598-018-24830-0 role in the innate immunity's specificity 24 . In initial cancer stages, CPA6 mRNA is more highly expressed than in late stages 25 . These genes are of importance to those looking for mechanisms of disease resistance or the genetic impact on vaccination efficacy.

Conclusions
To the authors' knowledge, this is the first time the Harderian gland transcriptome has been analyzed using RNA-seq. The Harderian gland of two genetic lines showed very different responses to NDV in number of DEG and amount of detectable NDV. The pathways activated in the nonchallenged Fayoumi suggest a more elevated immune system relative to the Leghorn under resting conditions. The challenged Fayoumi had significantly higher viral transcript counts in the Harderian gland at 2 dpi but no longer had detectable viral transcript counts at 6 dpi, unlike the Leghorn. This suggested the Fayoumis might clear the virus more quickly than the Leghorns. The Harderian gland is a first responder tissue after eye-drop vaccination or aerosol transmission of the virus. This study established a foundation for future research to investigate this tissue's unique role in host defense.

Animal experiment. This study was approved by the Iowa State University Institutional Animal Care and
Use Committee (IACUC log number 1-13-7490-G) and the experimental design and methodology has been described previously 17,18 . All experiments were performed in accordance with the committee's relevant guidelines and regulations. Two inbred lines, the Fayoumi (M 15.2) and Leghorn (GHs 6), from the Iowa State University Poultry Farm (Ames, IA), were used as a discovery platform to model relative resistance and susceptibility to NDV. Fayoumi and Leghorn breeders are a National Poultry Improvement Plan (NPIP) certified flock, tested and confirmed Salmonella and avian influenza virus free. Breeders of both lines are housed in the same building at the Iowa State University Poultry Farm (Ames, IA). At hatch, birds were placed in a biosafety level 2 facility. Chickens of both genetic lines were co-mingled throughout the study. Three weeks later, the challenged birds (n = 49) were inoculated with 200 μL of 10 7 EID 50 of La Sota NDV via the ocular/nasal route, and the nonchallenged birds (n = 40) were inoculated with PBS in the same manner. All challenged and nonchallenged birds were negative for NDV specific antibodies prior to challenge. Challenged and nonchallenged birds were kept in separate rooms from hatch until the end of the experiment. Approximately one third of the chickens were euthanized at 2, 6, and 10 dpi for tissue collection, the Harderian gland was removed, and placed into RNAlater (ThermoFisher Scientific, Waltham, MA). In the same week, tissues were processed and then stored at −80 °C. Each treatment group included four biological replicates for transcription analysis with the exception of the Leghorns at 2 dpi (5 challenged, 3 nonchallenged) and 3 challenged Fayoumis at 10 dpi. Total RNA was isolated from these samples using the RNAqueous kit (Thermo Fisher Scientific, Waltham, MA) and was DNase treated using a DNA-free kit (Thermo Fisher Scientific, Waltham, MA). All samples had an RNA quality number >8.0, as measured using the PROSize ® Data Analysis Software and the Fragment Analyzer ™ Automated CE System (Advanced Analytical Technologies, Inc., Ankeny, IA). The isolated RNA from each sample was used to independently construct two cDNA libraries for each biological replicate using the high-throughput protocol in the TruSeq RNA sample preparation guide (v2; Illumina, San Diego, CA). Libraries were validated using the Fragment Analyzer ™ Automated CE System (Advanced Analytical Technologies, Inc., Ankeny, IA) and 100 bp single end reads were sequenced on the HiSeq2500 (Illumina, San Diego, CA) at Iowa State University DNA Facility (Ames, IA). The sequencing data are available in the ArrayExpress database at EMBL-EBI (https://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-6038.
RNA-seq pipeline. The sequence data were taken through a standard pipeline for analysis using FASTX, TopHat2 26 , and HTSeq 27 as described previously 18 . A chi-square test was performed on the HTSeq output for each pair of technical replicates to confirm the technical replicates had minimal differences. Then the raw counts were combined for each pair of technical replicates. One technical replicate was discarded because only 10% of its reads mapped to the chicken reference genome, indicating technical failure.
Reads that did not map to the chicken reference genome were analyzed further as described previously 17 . These unmapped reads from each technical replicate were mapped to the La Sota genome (GenBank accession number JF950510.1) using BWA 28 and the mapped reads were counted with HTSeq 27 . The technical replicates' raw counts for each viral gene were summed for each biological replicate. The counts per million for each viral gene for each biological replicate were statistically analyzed in JMP statistical software (JMP Group Inc., San Francisco, CA) using a standard least squares, effect leverage full factorial test that accounted for line, dpi, gene, and their interactions.
Principal component analysis. For data visualization, pcaExplorer 29 was utilized. The dds function from DESeq2 30 was used for normalization, accounting for line, challenged, dpi, and sex in the model, and the rlogtransformation in pcaExplorer was used. The 5000 transcripts that contributed the most variance were used to calculate the variance associated with the principal components. Differential expression analysis. For differential expression analysis of the count data, the generalized linear model in edgeR 31 accounted for line, challenge, and dpi to determine the number of DEG (False Discovery Rate (FDR) <0.05). Contrasts were written to compare challenged and nonchallenged birds within each line and time, Fayoumis and Leghorns within each challenge group and time, and the challenge*line interaction at each time.
The edgeR output, including LFC and FDR for each transcript for each contrast, was input into Ingenuity Pathway Analysis (IPA; QIAGEN Inc., Redwood City, CA, https://www.qiagenbioinformatics.com/products/ ingenuity-pathway-analysis) 32 . The results from the challenged vs. nonchallenged Leghorns at 6 dpi contrast were input into IPA, transcripts with FDR <0.05 and absolute LFC >1 were used for an expression analysis. Figure 3 was generated based on the upstream regulator TCR (activation z-score = 2.712; p-value = 7.52E-10) from the upstream analysis portion in IPA. Additional diseases & functions (p-value <0.0001) relevant to NDV challenge were added to the network using the grow function (Fig. 3). The results from the Fayoumi vs. Leghorn contrasts were input into IPA and transcripts with FDR <0.10 were used for a comparison analysis for both challenged and nonchallenged birds at 2 and 6 dpi (Fig. 4); no pathways were significantly different between the two lines at 10 dpi. IPA calculates a z-score based on the expression levels of genes within a pathway or related to a functional term. If a z-score is positive IPA predicts activation of a pathway or function, and if negative IPA predicts inhibition of a pathway or function.

Method validation.
To validate the RNA-seq methodology, high-throughput qPCR using the Fluidigm Biomark system was performed. The methods 17,18 and primers 17 were reported previously. The LFC calculated by edgeR for the RNA-seq data and the −2 −ΔΔCT method from the Fluidigm Biomark output was used for the correlation of the challenged vs. nonchallenged comparisons at each time and within each line for 32 genes, resulting in 192 data points.