What Makes the Harderian Gland Transcriptome Different From Other Chicken Immune Tissues? A Gene Expression Comparative Analysis

The Harderian gland is a sparsely characterized immune tissue known to play an important role in local immunity. The function of the Harderian gland, however, is not clearly defined. Measuring the expression of all genes using RNA-seq enables the identification of genes, pathways, or networks of interest. Our relative RNA-seq expression analysis compared the chicken Harderian gland transcriptome to other important primary and secondary immune tissues including the bursa of Fabricius, thymus, and spleen of non-challenged birds. A total of 2,386 transcripts were identified as highly expressed in the Harderian gland. Gene set enrichment showed the importance of G-protein coupled receptor signaling and several immune pathways. Among the genes highly expressed in the Harderian gland were 48 miRNAs, a category of genetic elements involved in regulation of gene expression. Several identified miRNAs have immune related functions. This analysis gives insight to the unique immune processes inherent in the Harderian gland.


INTRODUCTION
Avian species have many unique immunological features compared to mammals with whom they last shared a common ancestor over 310 million years ago (Hedges, 2002). In birds the spleen is the largest lymphoid tissue, but is only able to encounter antigens that circulate through the blood because unlike mammals, birds lack a lymphatic system (Oláh et al., 2013). T cell development is similar in mammals and birds, but chickens have more γδ T cells than humans (Smith and Göbel, 2013). The thymus is a primary immune tissue where T cell development, differentiation, and maturation occurs.
Humans completely lack the bursa of Fabricius and only have a rudimentary Harderian gland, whereas these two tissues play very important roles in the chicken immune system. The bursa is a unique primary immune organ found in birds that plays a critical role in the immune response. B cell development, proliferation, and diversification occurs in the bursa, where B cells also undergo immunoglobulin rearrangement to create B cell receptors and mature B cells (Glick et al., 1956). The Harderian gland is located behind the eyes of the chicken and its function is not clearly defined, but includes the lubrication of the nictitating membrane (Bang and Bang, 1968). It is a relatively small tissue; in adult chickens the average weight was found to be 84.4 mg (Wight et al., 1971). The Harderian gland is known to contain many B cells. The majority of cells within the Harderian gland react to anti-chicken bursa cell serum (Albini and Wick, 1974). Lymphocytes from the bursa migrate to the Harderian gland prior to hatch and may not be involved in systemic immunity (Mueller et al., 1971;Baba et al., 1988). Also, terminal B cell maturation may occur in the Harderian gland (Manisikka et al., 1989). The Harderian gland is also home to large numbers of T cells. Equal numbers of CD3+ and Bu-1+ cells were found in the Harderian gland of both control and vaccinated chicks, and there were twice as many CD4+ than CD8+ cells in unvaccinated chicks as measured by immunostaining using monoclonal antibodies (Russell et al., 1997). The B and T cells of the Harderian gland play an important role in local immunity (Manisikka et al., 1989;Maslak, 1994). The bursa, thymus, spleen, and Harderian gland are among the most important immune tissues in the chicken.
Previously, these four immune tissues were compared directly via immunohistochemistry staining. In ducks, induction of CD8+ cell immunity within the spleen and thymus was stronger after challenge with an attenuated strain of hepatitis A, whereas in the bursa and Harderian gland CD8+ cell immunity was induced more strongly after challenge with the virulent strain (Ou et al., 2017). These tissues respond differently to antigen. In unstimulated chickens, µ Heavy chain and λ Light chain mRNA were expressed higher in the Harderian gland than the bursa, spleen, and thymus (Manisikka et al., 1989). Studying the transcriptome of these tissues elucidates the mechanisms utilized in response to pathogens. Until recently, the Harderian gland transcriptome had never been analyzed .
Unlike the Harderian gland, the transcriptomes of the bursa, spleen, and thymus tissues are well-characterized. Transcriptome analysis of the bursa revealed BCR receptor signaling and cytokine-cytokine receptor interaction pathways were impacted by avian pathogenic E. coli (APEC) (Sun et al., 2015), apoptosis of IgM+ cells, infiltration of macrophages, and increased expression of pro-inflammatory genes were seen after infection with velogenic Newcastle disease virus (NDV) (Kristeen-Teo et al., 2017), and defense response to virus, positive regulation of T cell-mediated cytotoxicity, and IFN-γ production pathways were impacted by infectious bursal disease virus (IBDV) infection (Ou et al., 2017).
The spleen transcriptome responded to a combined heat stress and LPS treatment by altering the expression of genes within the Hepatic Fibrosis/Hepatic Stellate Cell Activation and Macrophages pathway in two distinct genetic lines (Van Goor et al., 2017). In response to APEC infection, broiler splenic gene expression was predicted to affect the Jak-STAT and cytokinecytokine receptor signaling pathway (Sandford et al., 2011), and the spleen responded to NDV challenge by activating interferonstimulated genes (Zhang et al., 2018).
Compared to the bursa and spleen, there have been relatively few RNA-seq studies conducted on the chicken thymus. The thymus transcriptome responded to APEC by impacting the TLR signaling pathway, lysosome pathway, CAMs, and TCR signaling pathway (Sun et al., 2016). Another study showed thymus atrophy and its possible relationship with the expression of immune genes after challenge with LPS and Salmonella (Huang et al., 2016). In response to heat stress and an LPS challenge in the thymus transcriptome, ILK Signaling, Integrin Signaling, and cell proliferation pathways were all impacted (Lamont, personal communication).
Within each immune tissue, pathogen, strain, dose, time, genetic line, and more, greatly impact gene expression. Under basal conditions it is unclear how these tissues' transcriptomes compare, especially how they compare to the Harderian gland. A relative expression analysis of these fundamental immune tissues will help to better characterize the Harderian gland by identifying genes highly expressed (relative expression value greater than 2 SD from the mean) in this tissue relative to the bursa, thymus, and spleen. We assume the genes highly expressed in the Harderian gland are either related to tissuespecific non-immune function of the gland, or related to the unique immune function of this tissue in contrast to the other immune tissues studied. We hypothesize that the Harderian gland has mechanisms of defense that can be triggered rapidly because of its role in local immunity compared to the other more systemic immune tissues, and that the functional analysis of the genes highly expressed in the Harderian gland compared to the bursa, thymus, and spleen may elucidate these mechanisms.

Sample Descriptions and Processing
The Fayoumis (Line M 15.2) from the Iowa State University Poultry Farm (Ames, IA, United States) have been maintained as an inbred line since 1954 resulting in an inbreeding coefficient of 99.95% (Fleming et al., 2016). All publically available RNAseq data comes from the non-challenged Fayoumi controls from either a NDV challenge experiment (Experiment 1) Zhang et al., 2018) or a combined heat stress and LPS  Each bar represents the number of transcripts with a relative expression value within that range (x-axis). The red vertical lines represent two standard deviations below the mean (left) and above the mean (right). Transcripts to the left of the left red line are highly expressed in the other immune organs (spleen, bursa, and thymus). Transcripts to the right of the right red line are highly expressed in the Harderian gland. Figure Table 1).
In both experiments the Fayoumis were raised in floor pens with wood chips and ad libitum access to food and water. Although performed in separate batches, all tissues were collected and placed into RNAlater solution (Thermo Fisher Scientific, Waltham, MA, United States) for short-term storage, tissues were homogenized using mechanical disruption, RNA was isolated using an RNAqueous kit (Thermo Fisher Scientific, Waltham, MA, United States), DNAse treated with the DNA-free kit (Thermo Fisher Scientific, Waltham, MA, United States), and assessed for quality (RQN or RIN > 8). All samples underwent the same protocols to generate the cDNA libraries (TruSeq RNA sample preparation guide (v2; Illumina, San Diego, CA, United States), and were sequenced on the same HiSeq2500 machine to generate 100 bp single-end reads at the Iowa State University DNA Facility (Ames, IA, United States) ( Table 1). Spleen samples were collected in both Experiments 1 and 2. From Experiment 1 at ages 23 and 27 days, three of the four spleen and Harderian gland tissue samples were from the same individuals. No spleen samples from Experiment 1 were analyzed at 31 days of age. The spleen, thymus, and bursa samples from Experiment 2 came from the same four individuals ( Table 1).
The RNA-seq data underwent a standard pipeline previously described (Deist et al., 2017) and was mapped to the Gallus_gallus-5.0 (GCA_000002315.3) reference genome using TopHat2 (Kim et al., 2013). The number of reads mapped to each transcript was counted using HTSeq (Anders et al., 2015). All transcripts with less than four counts across all samples were removed resulting in 18,123 of 38,118 usable transcripts.

Calculating Relative Expression Values
The protocol for calculating relative expression values was previously described (Bailey et al., 2009;Pritchett et al., 2017). Pritchett et al. (2017) used the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) normalization method to normalize their counts. The current study has adapted the method using count data normalized with the variance stabilizing transformation in DESeq2 (fittype = mean; blind = true) (Love et al., 2014). A constant (2.32) was added to all normalized counts to make all values positive. The following formula was used to calculate the relative expression values (rEx). A transcript's rEx value was considered significant if it was more than two standard deviations from the mean. The same significance threshold was used previously (Pritchett et al., 2017). Comparing the maximum value to the median emphasizes the identification of transcripts highly expressed in the Harderian gland. Although this method may be sensitive to outliers, the standard deviation of individuals' normalized counts within each transcript in the Harderian gland was on average 1.05 (maximum SD = 5.07, minimum = 2. 26E-19). The individual sample from which the maximum normalized count value was obtained was well-represented across all samples. Each Harderian gland sample contributed a maximum normalized count value for at least 1,324 transcripts and at most 3,859 transcripts.

Analysis of Relative Expression Data
For data visualization, pcaExplorer (Marini, 2016) PCA plots were generated using dds and vst normalization from DESeq2 (Love et al., 2014) accounting for tissue and individual. The top 1000 most variant transcripts were used to calculate the principal components.
Transcripts highly expressed in the Harderian gland were further analyzed using Panther (Mi et al., 2017), Ingenuity Pathway Analysis (IPA; Qiagen, Redwood City, CA, United States), and STRING (Szklarczyk et al., 2017). These transcripts were converted to their associated gene name using BioMart on Ensembl (version 89) and input into Panther. Panther recognized 757 of the 992 input genes for an overrepresentation test using the GO biological process complete annotation set and the Gallus gallus reference list with a Bonferroni correction for multiple testing. Ensembl transcript identifiers (IDs) and the relative expression values for the transcripts highly expressed in the Harderian gland were used as input to IPA. Of the 2,386 transcripts, 942 were mapped FIGURE 3 | Significant GO terms associated with genes highly expressed in the Harderian gland. Transcripts highly expressed in the Harderian gland were converted to their corresponding gene ID using Ensembl (766 total) and input into Panther for statistical overrepresentation analysis using GO biological process (757 genes recognized and used for the analysis). The fold enrichment (purple bar) was calculated as the number of observed input genes divided by the number of expected genes based on the number of genes in the chicken genome. The number of genes associated with each GO term is listed to the right of the GO term in parenthesis. The Bonferroni correction for multiple testing was used to adjust p-values. The -Log 10 (p-value) for each GO term is represented by orange markers.
(identified) by IPA and used for analysis. Several canonical pathways were identified as significant, and those pathways with p-values less than 0.05 and included more than five genes have been reported. Associated gene names of transcripts highly expressed in the Harderian gland were input into STRING and used to generate a network. A high confidence (0.700) was used, all unconnected nodes were removed, and MCL clustering (inflation parameter = 3) was performed. A total of 836 nodes and 669 edges were included.

Principal Component Analysis
Samples clustered very tightly by tissue (Figure 1). Bursa, spleen, and thymus samples from Experiment 2 came from the same four birds, and there was also overlap between the birds that contributed the Harderian gland and spleen samples from Experiment 1 (Table 1). However, no clustering by individual bird was observed. Also, samples did not cluster by age or experiment. The large first principal component (70.3%) separated the Harderian gland from the other three tissues, showing the distinctiveness of the Harderian gland transcriptome. The second principal component separated the bursa, thymus, and spleen (PC2 = 13.9%). The two primary immune tissues, bursa and thymus, clustered more closely than the spleen (Figure 1). The spleen samples clustered very tightly together within group even though they were different ages and from different experiments (Figure 1 and Table 1). Six of the ten top loading genes for principal component 1 were also identified as highly expressed in the Harderian gland and included: MYL1, MYL2, CKMT2, and TMEM182. A list of the top and bottom loading genes was included in Supplementary Data Sheet S1.

Relative Expression Analysis
The relative expression values followed a tri-modal distribution (mean = 1.676; SD = 3.997) (Figure 2). Gaps in the distribution may be related to the distinct clustering of the Harderian gland and other immune tissues (PC1 = 70.3%; Figure 1). A total of 143 transcripts were highly expressed in the non-Harderian immune tissues, whereas 2,386 transcripts were highly expressed in the Harderian gland. The 2,386 transcripts were input into Ingenuity Pathway Analysis (IPA; Qiagen, Redwood City, CA, United States) and of the 936 identified by IPA, 96 were classified as transcription regulators and 4 as translation regulators. Of the 143 transcripts highly expressed in the other immune tissues only 53 had an associated gene name. No significant GO terms were identified. Some immune genes of interest in the 143 transcripts included C1QL3, C8A, and TLX1. For a complete list of the 143 transcripts see Supplementary Data Sheet S2. The relative expression analysis was more stringent than a differential expression analysis in which 99% of the transcripts were differentially expressed (false discovery rate < 0.05; data not shown). Gene Set Enrichment Analyses GO term analysis, pathway analysis, and network analysis were applied to the genes highly expressed in the Harderian gland on the assumption that these genes are the main drivers of functions that differentiate the Harderian gland from the other immune tissues. A cell type enrichment analysis (Shoemaker et al., 2012) showed no significant enrichment of any cell type based on these genes (data not shown); therefore, differences in expression levels among these tissues were likely not due to large differences in cell-type composition. The top-layer GO terms identified by Panther for the genes highly expressed in the Harderian gland are shown (Figure 3). Most GO terms were related to development and morphogenesis. The most significant GO term was G-protein coupled receptor signaling (Figure 3). No classic immune related GO terms were identified, however, cell fate commitment, G-protein coupled receptor signaling, and cell-cell signaling may be immune related.
IPA identified pathways associated with genes highly expressed in the Harderian gland ( Table 2). IPA identified more immune related pathways than the Panther GO term analysis (Figure 3). Notably, G-protein coupled receptor signaling was represented in both analyses (Figure 3 and Table 2). Calcium Signaling, ILK signaling, CXCR4 signaling, Crosstalk between Dendritic Cells and Natural Killer Cells, and MIF Regulation of FIGURE 4 | Network of genes highly expressed in the Harderian gland. Network generated by STRING. Edge thickness represents the confidence or strength of data support. A high confidence (0.700) cut-off was used to generate the network and all unconnected nodes were removed. Node color was based on MCL clustering (inflation parameter = 3). Supplementary Figure S1 includes the full network where a total of 836 nodes and 669 edges were included. The high number of nodes and edges led to a large network. From that network four clusters of interest (A-D) were chosen to display.
Innate Immunity are all pathways that have a direct relationship with the immune response.
The genes highly expressed in the Harderian gland formed a network with significantly more interactions than expected (p = 6.6e-16) (Supplementary Figure S1). STRING identified two significant (FDR < 0.05) KEGG pathways associated with these genes including Neuroactive ligand-receptor interaction (04080) and Tight junction (04530). Specific clusters of interest from the large network included Wnt genes (Figure 4A), GABA genes ( Figure 4B), Heat shock proteins (Figure 4C), and G-protein coupled receptors (GPCR) (Figure 4D).
A total of 48 miRNAs were highly expressed in the Harderian gland (Table 3), representing 39% of the miRNAs identified in this analysis. As the library construction kit used in this study employs a poly-A-tail selection, the reads that mapped to miRNA were likely from pre-processed miRNA. The correlation between the abundance of precursor miRNA and mature miRNA is dependent on the tissue and miRNA (Lee et al., 2008). It is possible that trace amounts of mature miRNA remained after poly-A-tail selection, but because the same kit was used for all tissues analyzed in this study, we assume no bias amongst tissues. It is also possible that reads belonging to the miRNA target sequence in regulated genes were incorrectly mapped to the miRNA itself. However, since reads that mapped to multiple places in the genome were removed, this is less likely.

DISCUSSION
The immune related genes that were highly expressed in the Harderian gland are of particular interest because they show how the homeostatic state is different compared to the bursa, spleen, and thymus. GPCR related genes, functions, and pathways were represented in every functional analysis. Every GPCR is activated by a specific ligand, which results in downstream signaling events. Glutamate is the ligand for several GPCRs identified as highly expressed in the Harderian gland, and IPA identified Glutamate Receptor Signaling as an impacted pathway. Glutamate is an amino acid that functions as a neuroimmuno-transmitter (Ganor and Levite, 2012). Glutamate can impact the immune system in several ways, i.e., affecting T cell survival, calcium levels, and cytokine expression levels (Ganor and Levite, 2012). Also, T cells, B cells, macrophages, and dendritic cells, express high levels of glutamate receptors (Ganor and Levite, 2012). The neurotransmitter GABA also acts as an immunomodulator (Jin et al., 2013). The genes highly expressed in the Harderian gland significantly impacted the GABA Receptor Signaling pathway and GABA receptors were represented in the network generated by STRING. GABA signaling can impact chemotaxis, phagocytosis, and cytokine secretion in immune cells (Jin et al., 2013). It is known there is a strong link between the central nervous system and the immune system (Black, 1994). The utilization of these neurotransmitters and their receptors in the Harderian gland compared to other immune tissues may be due to the focal role of local immunity in the Harderian gland, whereas, the other immune tissues in the current study are generally more important to systemic immunity. Another GPCR, FZD1, is the receptor for Wnt proteins. Several Wnt genes were identified as highly expressed in the Harderian gland and were represented in the network analysis. The Wnt pathway is tightly regulated, as it is involved in development, cell differentiation, and the immune response (Staal et al., 2008). A combined heat and LPS stress event showed the bursa transcriptome decreased expression of Wnt signaling genes (Lamont, personal communication). The miRNA mir-301a, has been shown to be activated by the Wnt pathway in glioma cells (Yue et al., 2016). Let-7 is also involved in cell proliferation and differentiation (Roush and Slack, 2008). Previously, increased levels of let-7 pri-miRNA were found in undifferentiated embryonic stem cells (Wulczyn et al., 2007). Increased expression of Wnt genes and let-7 miRNAs in the Harderian gland may increase the efficiency of this lymphoid tissue, by serving as a home to progenitor or naïve lymphocytes that can quickly be differentiated or activated in response to a stimulus.
Small, non-coding, miRNAs regulate gene expression levels post-transcription. Many of the miRNAs identified as highly expressed in the Harderian gland influence immune related genes and pathways. In human cell lines, mir-200b was shown to inhibit the TLR4 pathway (Wendlandt et al., 2012), mir-221 increased the expression of NF-κβ and STAT3 (Liu et al., 2014), and mir-193b overexpression resulted in increased autophagy (Nyhan et al., 2016). In chickens, mir-1764 decreased the expression levels of the inflammatory cytokine STAT1 (Jeong et al., 2013), and mir-30d may regulate IRF4 (Li et al., 2017). Avian influenza impacted the expression levels of several miRNAs in the chicken trachea and lung including the following miRNAs significant in the current study: let-7a-1, let-7f, let-7j, mir-1b, mir-30d, mir-34c, mir-101-2, mir-144, mir-200b, and mir-451 (Wang et al., 2009). Increased expression of these miRNAs in the Harderian gland would clearly impact the immune response in this tissue. Using miRNAs for regulation purposes in a tissue that is required to respond quickly to pathogens that enter via the eye is a useful strategy. Since the Harderian gland transcriptome has not been analyzed until recently  and was not used to generate the reference genome, there were likely several Harderian gland specific lncRNAs and miRNAs not identified in this analysis. Further research is necessary to confirm the correlation of abundance between the likely precursor miRNAs identified in this study and the mature miRNA levels, and to identify tissue specific, novel RNAs.
Overall this study identified 2,386 transcripts that were highly expressed in the Harderian gland compared to the bursa, thymus, and spleen of non-challenged chickens. These transcripts highlighted the interaction between the central nervous system and the immune system via the neuroimmuno-transmitters glutamate and GABA. The Harderian gland may utilize Wnt genes and let-7 miRNAs to control or stall cell differentiation. Previous studies have shown, several miRNAs identified as highly expressed in the Harderian gland have immune function. Moreover, these results elucidated the unique immune properties of the Harderian gland, a local immune tissue that must quickly respond to pathogens and vaccines that enter via the eye. It is important to gain a better understanding of the physiology of important avian tissues to develop better vaccines, breed for disease resistance in chickens, and potentially prevent zoonotic outbreaks.

ETHICS STATEMENT
The data used in this study was derived from animal subjects, but has been downloaded from public repositories. The ethics statements can be found in the manuscripts describing the original studies (Van Goor et al., 2017;Zhang et al., 2018).

AUTHOR CONTRIBUTIONS
MD: collected data, processed RNA-seq data, executed analysis, data interpretation, and wrote the manuscript. SL: oversaw data collection, analysis, interpretation, read and reviewed manuscript.

FUNDING
This study was funded by USAID Feed the Future Innovation Lab for Genomics to Improve Poultry and Hatch project 5357. MD was partially supported by the USDA NIFA 2013-38420-20496.

ACKNOWLEDGMENTS
The authors thank those who generated some of the data including Angelica Van Goor, Jibin Zhang, Melissa Monson, and Michael Kaiser as well as the entire Lamont lab group for their help processing tissues.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.00492/full#supplementary-material FIGURE S1 | The full network from which the clusters in Figure 4 were derived.
DATA SHEET S1 | A list of the top and bottom loading genes.
DATA SHEET S2 | For a complete list of the transcripts highly or lowly expressed in the Harderian gland compared to the bursa, thymus, and spleen.