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 (Deist et al., 2018).

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 cytokine–cytokine receptor signaling pathway (Sandford et al., 2011), and the spleen responded to NDV challenge by activating interferon-stimulated 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 tissue-specific 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.

Materials and Methods

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 RNA-seq data comes from the non-challenged Fayoumi controls from either a NDV challenge experiment (Experiment 1) (Deist et al., 2018; Zhang et al., 2018) or a combined heat stress and LPS experiment (Experiment 2) (Van Goor et al., 2017) (Table ). 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 ). 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 ).

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 (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.

Results

Principal Component Analysis

Samples clustered very tightly by tissue (Figure ). 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 ). 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 ). The spleen samples clustered very tightly together within group even though they were different ages and from different experiments (Figure and Table ). 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 .

Principal component analysis shows clear clustering by tissue type. Each dot represents a tissue from an individual chicken. Principal component 1 (PC1) separated the Harderian gland samples (green) from the other immune tissues. Principal component 2 (PC2) separated the Bursa (pink), Thymus (purple), and Spleen (blue). Ellipses were drawn with a 95% confidence. Dots with a black center were from Experiment 2. This plot was generated using pcaExplorer.

Relative Expression Analysis

The relative expression values followed a tri-modal distribution (mean = 1.676; SD = 3.997) (Figure ). Gaps in the distribution may be related to the distinct clustering of the Harderian gland and other immune tissues (PC1 = 70.3%; Figure ). 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 . 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).

Histogram of the relative expression values for each transcript. 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 generated in R.

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 ). Most GO terms were related to development and morphogenesis. The most significant GO term was G-protein coupled receptor signaling (Figure ). 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.

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 -Log10(p-value) for each GO term is represented by orange markers.

IPA identified pathways associated with genes highly expressed in the Harderian gland (Table ). IPA identified more immune related pathways than the Panther GO term analysis (Figure ). Notably, G-protein coupled receptor signaling was represented in both analyses (Figure and Table ). Calcium Signaling, ILK signaling, CXCR4 signaling, Crosstalk between Dendritic Cells and Natural Killer Cells, and MIF Regulation of 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 ). 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 ), GABA genes (Figure ), Heat shock proteins (Figure ), and G-protein coupled receptors (GPCR) (Figure ).

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 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.

A total of 48 miRNAs were highly expressed in the Harderian gland (Table ), 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 neuro-immuno-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 (Deist et al., 2018) 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 neuro-immuno-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.

Data Availability

The datasets analyzed in this study can be found at ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) or GEO (https://www.ncbi.nlm.nih.gov/geo) (see Table for more information).

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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer LB and handling Editor declared their shared affiliation.

Table 1

Sample information.

Tissue	Age (days)	Number of birds	Male:Female	Accessionc
Harderian glanda	23	4	2:2	E-MTAB-6038
Harderian glanda	27	4	3:1	E-MTAB-6038
Harderian glanda	31	4	2:2	E-MTAB-6038
Spleena	23	4	2:2	E-MTAB-5851
Spleena	27	4	3:1	E-MTAB-5851
Spleenb	22	4	1:3	GSE85434 (GEO)
Thymusb	22	4	1:3	E-MTAB-6290
Bursab	22	4	1:3	E-MTAB-6289

Table 2

Canonical pathways associated with genes highly expressed in the Harderian gland.

Ingenuity canonical pathways	p-Value	z-Score	Genes
Calcium signaling	0.0001	2.646	20
ILK signaling	0.0005	3.771	19
Thrombin signaling	0.0008	3	19
Actin cytoskeleton signaling	0.0028	3.873	19
Protein ubiquitination pathway	0.0141	—	19
Signaling by rho family GTPases	0.0072	4	18
G-Protein coupled receptor signaling	0.0347	—	18
Tight junction signaling	0.0006	—	17
cAMP-mediated signaling	0.0112	4	17
Cardiac hypertrophy signaling	0.0186	3.207	17
Phospholipase C signaling	0.0407	3.051	16
Cellular effects of sildenafil (Viagra)	0.0003	—	15
Epithelial adherens junction signaling	0.0010	—	15
RhoGDI signaling	0.0054	-3.464	15
Agranulocyte adhesion and diapedesis	0.0129	—	15
Mitochondrial dysfunction	0.0115	—	14
Sertoli cell-Sertoli cell junction signaling	0.0158	—	14
Hepatic fibrosis/hepatic stellate cell activation	0.0195	—	14
CREB signaling in neurons	0.0214	3.162	14
CXCR4 signaling	0.0195	1.897	13
Aldosterone signaling in epithelial cells	0.0219	—	13
GABA receptor signaling	0.0001	—	11
Gαi signaling	0.0107	1.897	11
Ovarian cancer signaling	0.0363	—	11
Synaptic long term depression	0.0407	3.317	11
Corticotropin releasing hormone signaling	0.0166	3	10
RhoA signaling	0.0324	3.162	9
Oxidative phosphorylation	0.0363	—	9
Transcriptional regulatory network in embryonic stem cells	0.0002	—	8
GPCR-mediated integration of enteroendocrine signaling exemplified…	0.0102	—	8
Crosstalk between dendritic cells and natural killer cells	0.0302	—	8
Regulation of actin-based motility by rho	0.0324	2.828	8
Glutamate receptor signaling	0.0087	—	7
Agrin interactions at neuromuscular junction	0.0234	2.449	7
Remodeling of epithelial adherens junctions	0.0234	—	7
Caveolar-mediated endocytosis signaling	0.0269	—	7
Basal cell carcinoma signaling	0.0288	1.633	7
Serotonin receptor signaling	0.0079	—	6
Ethanol degradation II	0.0174	—	5
Noradrenaline and adrenaline degradation	0.0234	—	5
MIF regulation of innate immunity	0.0309	2.236	5
Triacylglycerol biosynthesis	0.0398	—	5

Table 3

miRNAs highly expressed in the Harderian gland.

miRNA	Accessiona
gga-mir-1b	MI0001254
gga-let-7f	MI0001233
gga-let-7j	MI0001262
gga-let-7a-1	MI0001234
gga-mir-29b-2	MI0001266
gga-mir-29c	MI0001265
gga-mir-30c-2	MI0001205
gga-mir-30d	MI0001198
gga-mir-34c	MI0001261
gga-mir-101-2	MI0007558
gga-mir-130a	MI0001241
gga-mir-130b	MI0001239
gga-mir-133a-2	MI0001248
gga-mir-133b	MI0001206
gga-mir-133c	MI0001255
gga-mir-138-2	MI0001228
gga-mir-144	MI0004996
gga-mir-193b	MI0003698
gga-mir-200b	MI0001250
gga-mir-214	MI0008208
gga-mir-221	MI0001178
gga-mir-301a	MI0001240
gga-mir-365b	MI0022403
gga-mir-365-2	MI0003704
gga-mir-451	MI0004995
gga-mir-454	MI0006984
gga-mir-1592	MI0007319
gga-mir-1640	MI0007372
gga-mir-1684b	MI0022501
gga-mir-1737	MI0007476
gga-mir-1764	MI0007506
gga-mir-1772	MI0007515
gga-mir-1773	MI0007516
gga-mir-1800	MI0007544
gga-mir-1812	MI0007557
gga-mir-2126	MI0010731
gga-mir-6546	MI0022362
gga-mir-6580	MI0022399
gga-mir-6606	MI0022425
gga-mir-6609	MI0022428
gga-mir-6614	MI0022433
gga-mir-6653	MI0022473
gga-mir-6663	MI0022483
gga-mir-6668	MI0022488
gga-mir-6677	MI0022497
gga-mir-6701	MI0022523
gga-mir-6704	MI0022526
gga-mir-7474	MI0024147