The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by Enriching the Differentially Expressed Genes

The bursa of Fabricius (BF) is the central humoral immune organ unique to birds. The present study investigated the possible difference on a molecular level between two duck breeds. The digital gene expression profiling (DGE) technology was used to enrich the differentially expressed genes (DEGs) in BF between the Jianchang and Nonghua-P strains of ducks. DGE data identified 195 DEGs in the bursa. Gene Ontology (GO) analysis suggested that DEGs were mainly enriched in the metabolic pathways and ribosome components. Pathways analysis identified the spliceosome, RNA transport, RNA degradation process, Jak-STAT signaling pathway, TNF signaling pathway and B cell receptor signaling pathway. The results indicated that the main difference in the BF between the two duck strains was in the capabilities of protein formation and B cell development. These data have revealed the main divergence in the BF on a molecular level between genetically different duck breeds and may help to perform molecular breeding programs in poultry in the future.


INTRODUCTION
The bursa of Fabricius (BF) is the central humoral immune organ unique to birds, and it distinguishes the avian immune system from that of other species.In 1956, Bruce Glick revealed the importance of the BF in antibody production for the first time (Glick et al., 1956).It has been established that the BF plays a central role in B cell development and antibody production; in contrast to mammals where the immature B cells are mainly formed in the bone marrow (Osmond, 1986).The mature B lymphocyte can synthesize and secrete specific antibodies to participate in the humoral immune response when the body is under the stimulation of antigens (LeBien & Tedder, 2008;McHeyzer-Williams & McHeyzer-Williams, 2005).
Normal B cell development results from the interactions between the B cells and the supporting stromal cells in the BF (Nagy et al., 2004).Understanding B cell development in the context of the microenvironment is necessary to understand B cell biology.McCarthy et al. considered the entire BF organ as a model to investigate B cell biology in chickens using proteomics.This study identified 1753 B cell specific proteins, 1972 stroma specific proteins and 1473 proteins expressed in both tissue types (McCarthy et al., 2006).
Upon infection, mRNA expression of immune-related genes in the BF change significantly (Chen et al., 2013;Hang et al., 2014), suggesting that enrichment and clarification of immune-related genes in the BF, using a transcriptome technique, could help to reveal the functions of BF during avian resistance processes.Additionally, disease resistance breeding in poultry has become an important research field

The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by Enriching the Differentially Expressed Genes
in recent years (Thornton, 2010).Identifying candidate genes for disease resistance from different breeds will assist with breeding new disease resistant breeds.Next generation sequencing technology has provided new opportunities for helping to achieve these desired data and results The objective of this study was to identify the differentially expressed genes (DEGs) in the BF between two duck breeds, using the Digital gene expression profile (DGE) technology.One breed is a local species produced specially within the Sichuan Province in China, whereas the other breed is a newly cultivated breed from a highly artificial breeding process.Our results describe the fundamental processes involved in BF mediated disease resistance in ducks, providing key molecular information required to support future breeding programs in poultry.

Ethics Statement
All ducks were obtained from the Sichuan Agriculture University Waterfowl Breeding Experimental Farm, Sichuan, China.All of the procedures in this study were approved by the Animal Ethics Committee of Sichuan Agricultural University.All surgery was performed under sodium pentobarbital anesthesia, and all efforts were made to minimize suffering.

Duck Rearing and Sampling
Two genetically different duck breeds, the Jianchang (designated as JCH) and Nonghua P (NH-P) strains, were provided by Sichuan Agricultural University Waterfowl Breeding Experimental Farm.The ducks were hatched and maintained under identical conditions ad libitum until week 14. 6 embryos of duck were sacrificed at each stages of embryonic day 15, 20, and 25, and 6 ducks, half male and half female, were sacrificed every two weeks post hatching stages.The entire BF organ was weighed, and the BF organ index was calculated by the percentage of the BF weight compared to the corresponding bodyweights.Before being sacrificed, the ducks were starved for 12 h from the last meal, and then approximately 2 mL of blood was collected from the wing vein with EDTA (0.8 g/L) in a vacuum tube, and the plasma was then separated by centrifugation at 3000g for 10min at 4 °C.At weeks 3, 4, and 5, the BF organ tissues were isolated, and parts of them were immediately frozen in liquid nitrogen and then stored at -80 °C for RNA isolation.

ELISA
To show the innate immunity difference between JCH and NH-P ducks, the concentrations of IFN-γ, total protein and albumin were measured.The concentrations of IFN-γ, total protein and the albuminin blood serum were determined using ELISA kits (Qisong, Beijing, China) according to the manufacturer's instructions.Briefly, it was firstly to generate a standard curve after measuring the optical density of standard samples provided in the kits, then protein quantities of tested samples were calculated basing on theoptical density in the standard curve.The Globulin percentage was calculated by the following formula, Globulin percentage = (Total protein concentration -Albumin concentration)/ Total protein concentration*100%.

DGE Profiling Library Preparation and Sequencing
Two DGE libraries were constructed from the RNA pools extracted from BF tissues of JCH and NH-P ducks, respectively.Total RNA was extracted using Trizol reagent (Invitrogen, USA) according to the manufacturer's instructions.The quality of RNA was evaluated by the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA).The RNA integrity was assessed using the RNA Nano 6000 Assay Kit with the Bioanalyzer 2100 System (Agilent Technologies, CA, USA).Sequencing libraries were constructed using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, USA) according to the manufacturer's instructions, and index codes were added to attribute sequences to each sample.The library fragments were purified using the AMPure XP System (Beckman Coulter, Beverly, USA) to select cDNA fragments of preferentially 150-200 bp.PCR was performed using Phusion High-Fidelity DNA polymerase, universal PCR primers and an index (X) primer.The PCR products were purified (AMPure XP System), and the library quality was assessed using the Agilent Bioanalyzer 2100 System.Subsequently, the library preparations were sequenced on an Illumina Hiseq 2000 platform, and 100 bp single-end reads were generated.

Analysis and Mapping of DGE Reads
The sequenced data were filtered to remove raw reads containing adapter sequences, reads containing poly-N sequences and low-quality reads from raw data.In addition, the Q20, Q30 and GC content was calculated to evaluate the quality of data.The high quality reads were mapped to duck genome (BGI_ duck_1.0reference Annotation Release 101) database, allowing no more than 2 nucleotide mismatches.The high quality reads were designated as unambiguous

The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by
Enriching the Differentially Expressed Genes clean reads.For gene expression analysis, the number of unambiguous clean reads for each gene was calculated and normalized to RPKM (reads per kilobase transcriptome per million mapped reads).Prior to differential gene expression analysis, the read counts for each sequenced library were adjusted using the edge R software package through a one-scaling normalized factor.The differential expression analysis was performed using the DEGSeq R package (1.12.0) (Wang, Feng & Wang et al., 2010).A corrected P-value of 0.05 and a log2 (fold-change) of 0.05 were set as the threshold for significantly different expression.The Gene Ontology (GO) enrichment analysis of DEGs was implemented using the GOseq R package, in which the gene length bias was corrected (Young, Wakefield & Smyth et al., 2010).GO terms with corrected P-values less than 0.05 were considered to be significantly enriched with DEGs.KEGG is a database resource for understanding the high-level functions and utilities of biological systems (http://www.genome.jp/kegg/).The statistical enrichment of DEGs in KEGG pathways was examined using KOBAS software.

Validation of DEGs using real time PCR
Total RNA was extracted in BF from JCH and NH-P ducks.First-strand cDNA was obtained from the total RNA using a cDNA synthesis kit (Takara, Dalian, China) according to the manufacturer's instructions.The genes of interest were selected based on the list of DEGs identified from the DGE data (Table 3), and their sequence was downloaded from the published duck genome (BGI_duck_1.0reference Annotation Release 101).The GAPDH and β-actin was used as internal control genes.The primers were designed using Primer 5 software (Primer Biosoft International, USA), and the sequences are listed in Table 2.The mRNA expression levels of each gene of interest and the internal control gene were measured by real time PCR using the 96-well iCyclerIQ5 (USA) and the TaKaRa Ex Taq RT-PCR kit (Takara, Dalian, China).The following protocol was used: pre-denaturation at 95°C for 10 s, followed by 45 cycles at 95°C for 5 s and 62°C for 30 s.The annealing temperature and product length were different for each gene and are presented in Table 2.All reactions were performed in triplicate.The relative mRNA expression levels of the genes of interest were calculated using the "normalized relative quantification" method, followed by 2−ΔΔCT (Livak & Schmittgen, 2001).

Statistical analysis
All data were analyzed using the M.S. Excel program, and all results are listed as the mean ± S.D. A t-test was employed to analyze the significance of the differences between the groups at a level of p<0.05 using SPSS 13.0 software (Saddle River, NJ, USA).

Phenotypic divergence of BF between JCH and NH-P ducks
The BF organ index, serum concentration of IFN-γ and percentage of IgG in the total protein were investigated to detect the difference in immunity function of the BF between JCH and NH-P ducks.The BF organ index showed a clear divergence during the entire growth stage of the ducks.At the early embryonic stages (E15 and E20), the BF index in the NH-P ducks was slightly higher than in the JCH ducks, but at the later embryonic stages and post hatching stages, the BF index in the NH-P ducks was lower than in the JCH ducks (Figure 1A). Figure 1B&C shows the

The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by
Enriching the Differentially Expressed Genes differences in concentration of serum IFN-γ and IgG percentage in the total protein between the two duck breeds.Between the two types of ducks, the serum contents of IFN-γ showed no clear differences whereas the percentage of IgG in the total protein showed clear differences.A higher ratio of IgG was observed in the JCH ducks at early post hatching stages (W2-W8) and in the NH-P ducks at later post hatching stages (W10-W16).

DGE data quality and mapping reads to the reference database
Two DGE libraries (JCH and NH-P) were constructed and sequenced.There were approximately 18.9 and 14.3 million raw reads generated in the JCH and NH-P libraries, respectively (Table 1).More than 97% of the raw reads were filtered as clean reads, which were characterized as not containing N, were not low quantity and did not contain adapter sequence.Q20 and Q30 of the clean data were calculated, and the results showed that these items were more than 90% and 96%, respectively for the JCH and NH-P groups.The clean reads mapped to genes with no more than two nucleotide mismatches were near 12.7 and 9.7 million in the two individual libraries, more than 67% of which are perfectly matched (Table 2).
Saturation analysis was performed to confirm whether the number of detected genes increases proportionally to a sequencing coverage (total clean reads number).The data indicated that the DEGs with a RPKM value of more than 3-15 could be easily enriched using in-depth sequencing, and only small parts of genes with a RPKM value of less than 3-15could not be wholly enriched (Figure 2A & B).The uniform distribution suggested that the reads did not preferentially distribute near the 5' and 3' end of transcripts in the highly expressed transcripts whereas the reads showed a preferential distribution in the medium and lowly expressed transcripts (Figure 2C & D).

Analysis of differential gene expression and annotation
To enrich the genes with a significant difference in expression between JCH and NH-P, we calculated the gene expression levels using the fold changes with the default value calculated as |log2(FoldChange)|>0.05.The volcano plots shows the ranges of distribution in DEGs as a whole (Figure 3).A total of 95 up-regulated and 100 down-regulated genes were identified in JCH relative to NH-P with p-value ≤0.05All the DEGs enriched have been annotated (Table 2).The enrichment of DEGs in the BF were investigated to show the main possible differences in physiological processes between JCH and NH-P ducks.All of the enriched GO terms were presented in Table 3 (Tables  3, 4 and 5 are accessed by a direct contact with the author).Figure 4 A listed all significantly enriched GO terms.Most of the terms were distributed in the classifications of the cellular component, and fewer were distributed in the classifications of molecular function.Figure 4B showed the number of up and down-regulated genes in the significantly enriched GO terms in JCH ducks compared to NH-P ducks.The gene expression was up regulated in the majority of the GO terms of cellular components, in half of the GO terms of biological function, and in all of the GO terms of molecular function.Note: a The number of all clean reads ratio in total raw reads; b The number of all mapped reads ratio in total clean reads.Q20 represent the error rate less than 0.01, Q30 represent the error rate less than 0.001.

KEGG pathway enrichment in DEGs
Figure 5 lists the top 20 highly enriched pathways, and the results showed that the ribosome-related pathways have the maximum enrichment factor value and the lowest Q-value.Other pathways, such as spliceosome-related pathways, RNA transport and the RNA degradation process, were all enriched significantly by the DEGs.In all of the enriched pathways (Table 5), some of the DEGs were enriched in disease-related pathways, such as human T-cell leukemia virus type I (HTLV-I) infection, the B cell receptor signaling pathway, the T cell receptor signaling pathway, bacterial invasion of epithelial cells, antigen processing and presentation, TNF signaling pathways, the intestinal immune network for IgA production and the NF-κB signaling pathway.Some of the DEGs were enriched in the MAPK signaling pathway, the estrogen signaling pathway, apoptosis, the Jak-STAT signaling pathway, focal adhesion, and ECM-receptor interaction.

The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by Enriching the Differentially Expressed Genes
Verification of the DEGs Ten genes were selected in chance in the DEGs to verify the expression profiles.The up regulated genes included ANKRD11, CD3D, and THYN1, and the downregulated genes included DOCK8, HTT, IL17RA, MSN, MYBL1, NFATC3, and XPO1.The expression profiles of these selected genes were determined by real time PCR in the BF of JCH and NH-P ducks at weeks 3, 4, and 5.The results showed that 7 of these genes were consistent with the RPKM results in sequencing data at week 3.There were only 3 genes that exhibited similar expression trends between the two types of ducks, at week 4 and week 5, respectively (Figure 6).

DISCUSSION
The BF is unique to birds and plays a central role in avian immunity resistance.Our study compared the developmental characteristics and immune differences of JCH and NH-P breeds, by observing the growth states of BF and determining the serum IFN-γ and IgG protein concentrations.These data confirmed that differences in innate immunity existed in the BF between the duck breeds.Most likely due to the advantages of placing emphasis on economic traits and high selection pressures, some cultivated breeds lose important recessive genes, which may decrease immune resistance more than in local ducks.Similar data have been observed in chickens.Alshwabkeh et al. demonstrated that the Balady chicks (local breed) were more resistant than three commercial strains bred in the US (Hubbard), Holland (Lohman), and Canada (Shaver), when each chick was orally infected with 10(6) CFU of Salmonella gallinarum (Alshwabkeh, 2001).
A variety of factors, such as viruses, nutrients, etc., could affect the physiology and development status of the BF and affect changes at the molecular level (Guo et al., 2012;Lu et al., 2010;Peng et al., 2011;Yin et al., 2015), which supports complex molecular networks that regulate BF function in poultry.The present study identified a total of 195 DEGs between the JCH and NH-P ducks.This is the first list of differences of gene expression in the BF between two genetically different breeds.However, these DEGs were confirmed by a relatively low repeat ability of real time PCR data in this study, which was 70% at week 3 and only 30% at weeks 4 and 5.These discrepancies may be attributed to complex factors, such as the differences between the techniques of DGE and real time PCR.In this study, by DGE, the intermediate data for each step was strictly followed using a quality control analysis (Table 1, Figure 2), and all of these analyses showed that the DGE data were stable and accurate.
The BF functions as a lymph gland and plays a key role in the differentiation of B-lymphocytes during the first two to three months after the chicken hatches (Glick et al., 1956;Schat & Skinner, 2014).Figure 1A showed the duck BF begins to atrophy at week 4, using the BF index (Cazaban et al., 2015).The maximum size of the BF in ducks was observed slightly earlier than in chicks.In chicks, it was reported that the BF development begins during incubation and reaches its maximum size between 8 to 10 weeks of age (Oláh & Vervelde, 2008).In the present study, the BF of ducks was at the earliest stage of the atrophy process, and it showed differences between breeds.In the top 20 significantly DEGs, some of them are involved in apoptosis, such as TP53BP2 (Gorina & Pavletich, 1996) and ACTG2, and some of them take part in the regulation processes of gene expression, translation and protein modification, such as Med13l, DPP7, NELFA, Ctbp2, SERPINH1 (Christiansen, Schwarze & Pyott et al., 2010), PLEKHG4 and HTT.Some of the genes are related to the biological activities of B cells, such as CHRNB4, HFE (Drakesmith et al., 2002), PLEKHG4 and CD3D (Linderoth et al., 2008).These enriched DEGs are involved in different types of biological processes, which primarily supported the notion that multiple differences may exist in the BF between JCH and NH-P.
GO analyses were further employed to show the differences in molecular function, biological processes, and cellular components in the BF between the duck breeds.In detail, the differences were mainly in

The Gene Ontology Differs in Bursa of Fabricius Between Two Breeds of Ducks Post Hatching by Enriching the Differentially Expressed Genes
ribosome-related processes, structural constituents of the ribosome and gene expression and translation.The BF is an important lymph gland where B cell differentiation occurs in birds.The B cells are responsible for matching and detecting the antigens present in the body, and then secreting antibodies (LeBien & Tedder, 2008).Therefore, it was reasonable that the processes of gene expression and translation were up regulated, which may contribute to increased protein synthesis.In particular, we also observed some B cell activityrelated GO terms, such as B cell activation involved in the immune response, B cell cytokine production, B cell proliferation and B cell cytokine production.We confirmed that the BF can provide a useful experimental model of the development of B lymphocytes (Kincade, 1987).Korte et al. analyzed the proteome dynamics of the chicken BF organ during embryonic and post hatch development, and they have found that genes such as retinal-binding protein 5, the actin-cytoskeleton, vinculin, gelsolin, and chloride intracellular channel 2 are most likely critical to biological processes such as the migration of B cell precursors into the bursa anlage (Korte et al., 2013).
The DEGs were enriched in the KEGG database.The results still revealed that the main divergences referred to processes of the ribosome, spliceosome, RNA transport and RNA degradation, suggesting that protein synthesis was possibly the main difference between the two genetically different duck breeds.Luna-Acosta et al. observed the anti-apoptotic effect of growth factor (GH) in the BF, and they proved that the effects might be mediated through the PI3K/Akt pathway (Luna-Acosta, Alba-Betancourt and Martinez-Moreno et al., 2015).Additionally, Cheng et al. demonstrated that the NFκB-iNOS-NO signaling pathway may be important in the BF, due to its involvement in the process of generating for sythiaside, which exerts an anti-inflammatory effect (Cheng, Zhao & Li et al., 2014).NF-κB is also involved in B cell development (Franzoso, Carlson & Xing et al., 1997).In the present study (Table 4; Tables 2, 3 and 4 are accessed by a direct contact with the author), both the PI3K/Akt (matched 5 DEGs) and NF-κB (matched 2 DEGs) pathways were enriched, although they exhibited a low enrichment factor and high p-value, indicating that these two pathways play essential roles in the BF and may contribute less to the differences in the immunity of the BF between the JCH and NH-P ducks.Other pathways such as human T-cell leukemia virus type I (HTLV-I) infection, the B cell receptor signaling pathway, the T cell receptor signaling pathway, bacterial invasion of epithelial cells, antigen processing and presentation, TNF signaling pathways, and the intestinal immune network for IgA production were enriched in the present study.These pathways are related to immunity resistance and processes (Davis et al., 2010;Matsuoka, 2005;Schaeffer, Debnath & Yap et al., 1999), suggesting that there are differences in the immune reaction and processes in the BF between the two types of genetically different ducks.Some enriched pathways are highly related to the cell cycle and activities, such as the MAPK signaling pathway (Zhang & Liu, 2002), the estrogen signaling pathway (Foster et al., 2001), apoptosis, the Jak-STAT signaling pathway (Steelman Pohnert & Shelton et al., 2004), focal adhesion (Zhao, 1998), and ECM-receptor interaction (Lukashev & Werb, 1998), supporting the notion that the differences in the BF may also occur in cell number and interactive effects among cells between JCH and NH-P ducks.

CONCLUSION
195 DEGs in the bursa have been screened through DGE technology between the Jianchang and Nonghua-P duck breeds.These highly matched GO terms were concentrated in ribosome-related processes, structural constituents of the ribosome and gene expression and translation, and the KEGG pathways were mainly distributed in the ribosome, spliceosome, RNA transport and RNA degradation.These findings suggest that the main difference in the BF between Jianchang and Nonghua-P duck strains is most likely in the capabilities of protein formation and B cell development.

Figure 1 -
Figure 1 -Phenotypic divergence of bursa of Fabricius (BF) between JCH and NH-P ducks.A, A comparison of the BF index during embryonic stages and post hatching stages.B, A comparison of serum IFN-γ post hatching stages.C, A comparison of serum IGG in the percentage of total proteins in post hatching stages.At least 6 individuals for each breed and stage were employed to generate the data.

Figure 2 -
Figure 2 -The evaluation of sequencing quality.A & B, The saturation curve analysis in JCH and NH-P; C & D, The destiny distribution of the transcripts at different expression levels in JCH and NH-P.

Figure 5 -
Figure5-KEGG scatter plot analysis of DEGs.A total of 195 DEGs were screened in the KEGG databases to enrich the top 20 pathways differentially expressed in the BF of JCH and NH-P ducks.The smaller the q value, indicating that the signaling pathway is more significant, the greater the enrichment factor, the more DEGs that were enriched.

Figure 4 -
Figure 4 -GO enrichment analysis.A, All significantly enriched GO terms with a corrected p-value less than 0.05 by the DEGs.BP, biological process; CC, cellular component; MF, molecular function; B, The number of genes up and down regulated in the subgraph of the enriched GO terms, which were compared in JCH ducks and NH-P ducks.

Figure 6 -
Figure 6 -Validation of DEGs by real time PCR

Table 1 -
Primers used for qRT-PCR

Table 2 -
Quality analyses of digital gene expression (DGE) Library