Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Transcriptome profile in bursa of Fabricius reveals potential mode for stress-influenced immune function in chicken stress model

BMC Genomics volume 19, Article number: 918 (2018) Cite this article

Abstract

Background

The molecular mechanisms underlying stress-influenced immune function of chicken (Gallus Gallus) are not clear. The stress models can be established effectively by feeding chickens corticosterone (CORT) hormone. The bursa of Fabricius is a unique central immune organ of birds. RNA-Seq technology was used to investigate differences in the expression profiles of immune-related genes and associated pathways in the bursa of Fabricius to clarify molecular mechanisms. The aim of this study was to broaden the understanding of the stress-influenced immune function in chickens.

Results

Differentially expressed genes (DEGs) in the bursa of Fabricius between experimental group (basal diet with added CORT 30 mg/kg; C_B group) and control group (basal diet; B_B group) were identified by using RNA-seq technology. In total, we found 1434 significant DEGs (SDEGs), which included 199 upregulated and 1235 downregulated genes in the C_B group compared with the B_B group. The immune system process GO term was the top significantly GO term, including MYD88, TLR4, IL15, VEGFA gene and so on. The cytokine-cytokine receptor interaction pathway and the Toll-like receptor signaling pathway were the key pathways affected by stress. The protein-protein interaction (PPI) analysis of the SDEGs showed that VEGFA, MyD88 and IL15 were hub genes and module analysis showed that MYD88, TLR4 and VEGFA play important roles in response to stress.

Conclusion

This study showed that the VEGFA and ILs (such as IL15) via the cytokine-cytokine receptor interaction pathway, MYD88 and TLR4 via the Toll-like receptor signaling pathway may play important roles in the regulation of immune function under stress condition with CORT administration. The results of this study provide a reference for further studies of the molecular mechanisms of stress-influenced immune function.

Background

Stress is a physiological manifestation of the body’s defence against adverse environmental effects. It is nonspecific responses of the body to adapt to the environment and maintain homeostatic equilibrium balance [1]. With the rapid development of the modern poultry industry and highly intensive production management, more and more stress factors have appeared in the production process, such as damp-heat, congestion, transportation, conventional vaccination, rapid cultivation of high production breeds and so on. Moderate stress can promote the immunity of the body, but excessive stress not only affects the growth, production and reproduction performance but also affects the immune system of animals, causes immune suppression, diseases and even death [2]. Stress can increase the level or secretion of adrenaline; increase the blood pressure; increase the level of corticosteroids in the plasma; contribute to the degeneration of the immune organs such as thymus, bursa of Fabricius, spleen and lymphatic tissue [3]; change the number of leukocytes [4], thereby affecting the presence of the eosinophilic and heterophilic leukocytes [5]; cause T lymphocyte injury and macrophage inhibition; and increase the catabolism of Immunoglobulin G [6], thereby resulting in a decrease in immune response and a high susceptibility to epidemics. The influence of stress on poultry production is becoming more and more serious, which has brought huge economic losses to the modern poultry industry. For a long time, researchers and poultry producers have done researches and practical works on the prevention and control of production stress [7,8,9,10], but they have not fundamentally solved this problem.

Many stress factors are present during the intensive production of poultry, and single stress factors are not representative. A stress model, therefore, is an important means for studying the stress problem. It has an irreplaceable role in the study of the mechanism of stress and the response to stress, the avoidance of stress, and in research on the drugs that interfere with stress. Corticosterone, which is secreted by the adrenal cortex and is the main glucocorticoid of birds, is involved in the stress response of the body; regulates the metabolism of carbohydrate [11], lipid [12] and protein [13, 14]; maintains the normal activities of various tissues and organs; has roles in anti-inflammatory and inhibition of immune function; and has reduced humoral and cellular immune functions [11, 15]. The establishment of the chicken stress model with exogenous CORT is theoretically based and is feasible [16].

Researches on the effect of various stress on chicken immune function in local or commercial breeds is mainly evaluated by the changes in relative lymphoid organ weight [17], partial immune related indices [18], cytokines levels [19], immune cell apoptosis [20, 21] and so on. RNA-Seq has been used to analysis the transcriptomes of the different cells and tissues from various animals [22,23,24]. Applying transcriptome technology to study specific stress were paid more attention to growth and development of chicken [25, 26]. Studies related to poultry stress affecting immune organ function using RNA-Seq technology has been not reported. The bursa of Fabricius is the unique central immune organ of birds and is the location of B lymphocyte differentiation and maturation. Fully differentiated B lymphocytes migrate from the bursa of Fabricius to peripheral lymphoid organs to colonize, reproduce and perform important immune functions. B lymphocytes can be converted into plasma cells and produce antibodies that participate in humoral immunity. B lymphocytes are the only cells in the body that produce antibodies that play a central role in the immune response process [27].

Therefore, to analysis molecular regulation mechanisms of the stress influencing immune function, we administered CORT in an experimental diet that included 30 mg CORT/kg basal diet to build a stress model according to previous studies [16]. Then, the bursas of Fabricius from the chicken models were analyzed via transcriptome profiling by RNA-seq technology to deeply explore genes or pathways involved in the stress response. Verifying the function of these key genes involved in the immunoregulation can lay a foundation for the analysis of molecular regulation mechanisms of the stress influencing immune function.

Results

Illumina sequencing and analysis of mRNA expression in the bursa of Fabricius

Six libraries B_B_1, B_B_2 and B_B_3 from the control group and C_B_1, C_B_2 and C_B_3 from the experimental group were sequenced on an Illumina HiSeq platform, with 343 million reads in total being generated, of which 95.97% (329 million) passed the filter for clean reads. Of these, 85.22 to 86.94% clean reads were aligned to the chicken genome (Gallus gallus 4.0), the multiple mapped reads ranged from 1.84 to 1.90%, and the uniquely mapped reads ranged from 83.32 to 85.04%. Among these mapped reads, 68.1 to 72.0% were mapped in annotated exons, 19.3 to 22.9% were mapped in the intergenic regions, and 8.3 to 9.0% were mapped in introns. Most of reads were mapped in the same coding regions as observed in the previous study [22]. The GC contents of the clean reads was 48.33 to 49.7% (Additional file 1: Table S1). To test the reliability of the samples and whether the sample selection was reasonable, we performed a correlation analysis of the 6 samples. The results (Additional file 2: Fig. S1) showed that the intra-group correlation coefficient R2 was approximately 0.99, which was higher than the correlation coefficient R2 of approximately 0.93 between the groups (R2 > 0.92 indicated that the experimentation was relatively reasonable).

Analysis of gene expression levels in the bursa of Fabricius and DEGs introduced by CORT added to feed

Compared the sequencing results of two groups with the reference gene group, 19,254 genes were identified, of which 18,684 and 18,814 genes were identified in the experimental group and control group, respectively. Further analysis of the transcriptome data revealed that there were 440 and 570 genes specifically expressed in the experimental and control groups, respectively (Additional file 3: Fig. S2). Among the 19,254 genes, 15,832 had annotated information, whereas 3422 were unknown genes. The most highly expressed gene in the C_B group was cathelicidin-B1-like (CATHB1), and the gene was more highly expressed than in the B_B group although not significant difference. The highest expressed genes in the B_B group were the β-2-microglobulin (B2M, padj = 1.21E-86, log2FC = − 1.08) and immunoglobulin lambda-like polypeptide 1 (IGLL1, padj = 1.21E-51, log2FC = − 1.18), which were associated with immune response and were significantly upregulated expressed than in the C_B group (Additional file 4: Table S2).

In total, we found 1434 SDEGs, which included 199 upregulated and 1235 downregulated genes in the C_B group compared with the B_B group (Fig. 1; Additional file 5: Table S3). Herein, we marked the top 20 up- and downregulated genes between the C_B and B_B group according to the log2FC in Additional file 5: Table S3. Significantly downregulated genes included the radical S-adenosyl methionine domain-containing 2 (RSAD2); chemokine (C-C motif) ligand 19 (CCL19); chemokine-like ligand 1 precursor (CCL4); immunoresponsive 1 homolog (IRG1); IL4I1 and CCLI10. Significantly upregulated genes included the heterogeneous nuclear ribonucleoprotein K (HNRPK); Gallinacin-2 (GAL2) and CACNB4. The heat map of cluster analysis of SDEGs showed that the gene expression patterns were similar intra-group, while were different between groups (Fig. 2).

Fig. 1
figure 1

Comparison of SDEGs between the two groups. Scatter plot showing the correlation of gene abundance. SDEGs are shown as red (up) or blue (down) dots and gray dots indicate genes that lack of significantly difference

Full size image
Fig. 2
figure 2

Cluster analysis of SDEGs by the FPKM value. The X-axis indicates the samples in the different groups. The three samples on the left are from the C_B group, and the three samples on the right are from the B_B group. The Y-axis is the gene cluster across the C_B and B_B group. Red indicates the high expression genes, and blue indicates the low expression genes by the value of log10(FPKM + 1)

Full size image

Gene ontology (GO) enrichment analysis

Among the 1434 SDEGs, 1269 SDEGs of the C_B group/B_B group were enriched for 1052 GO terms (Additional file 6: Table S4). The top 30 most significantly enriched GO terms are shown in Fig. 3 (For details see Additional file 6: Table S4), which were all biological process terms (including immune system process, immune response, regulation of immune system process and so on). The immune system process was the most significantly enriched GO term. Totally, 243 SDEGs were enriched in the immune system process, including VEGFA, ILs (such as IL-1 β, IL8, IL7, IL18), THBS1, TLR4, TLR3, MYD88, ITGB3 and other genes, which may play important roles in the process of stress affecting immune function.

Fig. 3
figure 3

The top 30 significantly enriched GO terms of the SDEGs. The X-axis indicates the number of DEGs for each GO terms; the y-axis corresponds to the GO terms

Full size image

The Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis

KEGG pathway enrichment was used to identify the significantly biochemical metabolic pathways and signal transduction pathways modulated by CORT. The SDEGs of the C_B group/B_B group were annotated into 127 KEGG pathways. The top 20 enriched pathways of the SDEGs for the C_B and B_B group are shown in Fig. 4 (For details see Additional file 7: Table S5). In total, 9 KEGG pathways were significantly enriched, including the cytokine-cytokine receptor interaction pathway, the cell adhesion molecules pathway, the phagosome pathway, the influenza A pathway, the Toll-like receptor signaling pathway, the ECM-receptor interaction pathway, the NOD-like receptor signaling pathway, the focal adhesion pathway and the salmonella infection pathway.

Fig. 4
figure 4

The top 20 significantly enriched KEGG pathway of the SDEGs. The x-axis shows the rich factor; the y-axis corresponds to the KEGG pathways. The dot color represents the q-value, and the dot size represents the number of SDEGs enriched in the reference pathway

Full size image

The cytokine-cytokine receptor interaction pathway was the most significantly enriched pathway. We found that 48 SDEGs were enriched in the cytokine-cytokine receptor interaction pathway for the C_B vs B_B comparisons, including vascular endothelial growth factor (VEGFA), IL8, IL18, IL7, IL-1 β, IL15, CCL4, IL1R2 and so on. The Pearson’s correlation coefficient was 0.92 between the expression of the receptor genes and ligand genes in the cytokine-cytokine receptor interaction pathway (Additional file 8: Table S6). Interestingly, we found that the 48 genes enriched in the cytokine-cytokine receptor interaction pathway were all downregulated. Similarly, the 21 genes (including myeloid differentiation primary response protein 88 (MYD88), IL-1 β, TLR4, k60, IL8 and so on) were enriched in the Toll-like receptor signaling pathway and all but one gene (TLR1LA) were also downregulated in C_B group.

The protein-protein interaction (PPI) network analysis

The PPI networks for SDEGs consisted of 250 proteins and 341 pairs of PPIs (Additional file 9: Table S7). The hub proteins exceeded 8 interactions, including VEGFA, MYD88, IL15, SDC4, ITGB4, ITGB3, ITGB5, THBS1, JUN, k60, LAPR1, PLCG1, MAPK11, RHOC, and RHOJ (Additional file 10: Fig. S3). Only two proteins exceeded 13 interactions, ITGB3 and VEGFA, of which the VEGFA has 18 interactions. A total of 21 modules including 81 proteins were obtained using default criteria, in which 19 modules with node > 3 and a MCODE score > 3 (Additional file 11: Fig. S4). Among them, we found module 19th had the MYD88, TLR4 and BIRC2 genes, of which the MYD88 and TLR4 were enriched in the Toll-like receptor signaling pathway. The VEGFA gene were in module second and enriched in the cytokine-cytokine receptor interaction pathway (Fig. 5).

Fig. 5
figure 5

Modules of the protein-protein interaction network. The color of nodes indicates the degree of foldchange, the value was “-log2FoldChange”. The edge between proteins represent interactions between them and the confidence score were more than 0.7

Full size image

Validation of gene expression by quantitative real-time PCR (qRT-PCR)

To confirm the accuracy of the differential expression analysis results, qRT-PCR was performed to validate the 12 randomly selected DEGs, including C7, POLDIP2, CASP3 and CD79B, from the upregulated genes and HSP70, IDUA, MYD88, TGM2, BCL2, PCBD2, HYAL1 and IGSF3 from the downregulated genes, with GAPDH and β-action as the reference gene respectively. The qRT-PCR samples were the same as the RNA-seq samples, with three biological repeats for each group and three repetitions for each sample. The qRT-PCR results, with GAPDH and β-action as the reference genes, were consistent with the RNA-seq results (Fig. 6).

Fig. 6
figure 6

Compared the gene expression levels of RNA-seq with qRT-PCR. The left axis represents the expression levels determined by RNA-seq in FPKM units, and the right axis represents gene expression levels determined by qRT-PCR. Bars represent the mean (±SE) of three samples. The black column indicates the FPKM value; the grey column indicates the qRT-PCR used GAPDH as reference gene; the white column indicates the the qRT-PCR used β-action as reference gene

Full size image

Discussion

Glucocorticoid is a major stress hormone which increase rapidly in response to stress, such as CORT in bird [28]. Glucocorticoids regulate their effects by binding and activating glucocorticoid receptor (GR). Glucocorticoid receptor is a ligand-activated transcription factor, and has a zinc-finger motif DNA-binding domain, which can be used to regulate transcription both positively (transactivation) and negatively (trans-repression) [29]. In transactivation, GR activates downstream gene transcription by binding to a glucocorticoid response element (GRE). In trans-repression, GR interacts with other transcription factors, such as activator protein-1 (AP-1) and nuclear factor-kappa B (NF-κB), and indirectly inhibits transcription by reverse inhibition [30, 31]. Many of the GRs of the target genes are involved in regulating immune system and inflammatory response. The trans-repression activity of GR is the major basis for the anti-inflammatory and immunosuppressive effects of glucocorticoid. Previous research reported that feeding an experimental diet with 30 mg of added CORT/kg basal diet increased the plasma concentrations of CORT in chicken, indicating successful induction of stress response [16]. Therefore, we established a stress model by adding CORT to the diet at the same dose to analyze the molecular mechanism of stress affecting immune function [16].

In this study, we construct the stress model of chicken with administration of adding 30 mg CORT per kilogram of the basal diets referring to previous studies [16]. The aim of this study was to analysis molecular regulation mechanisms of the stress influencing immune function. Hence, the stress model was evaluated by detecting the indicators related to stress, including the physiological indexes such as CORT, glucose (GLU) and total protein (TP) in peripheral blood of chickens, serum CD3+, CD4+, IgG, TNF-α, IL-1β, IL-6 and relative lymphoid organ weights. The results showed that the levels of serum CORT, GLU, TP, TNF-α, IL-1β and IL-6 were increased significantly, the levels of alkaline phosphatase (ALP), CD3+ and CD4+ were decreased significantly (Additional files 12 and 13: Tables S8 and S9), however the relative bursa of Fabricius weights were no significant difference (results not shown). Corticosterone is a major stress hormone which increase rapidly in response to stress [32]. The higher serum CORT level, the more obvious the stress. The increase of serum GLU level after stress may be related to the induction of CORT. Corticosterone could reduce the utilization of GLU in various tissues of the body and enhance the hepatic gluconeogenesis. Moreover, CORT could promote protein metabolism and increases blood protein concentration. Cytokines are important mediators between the neuroendocrine system and the immune system. Stress could stimulate the release of inflammatory factors such as TNF-α, IL-1β, IL-6, and cause inflammatory response. Stress also could decrease the number of T cells and reduce cellular immunity function [33]. Our results indicated that the stress model of chicken was successfully constructed.

The bursa of Fabricius is the unique central immune organ of birds and play a central role in the immune response process [27]. B lymphocytes of birds are derived from bone marrow pluripotent stem cells and migrated to bursa of Fabricius, then in which differentiated and maturated. There are not only B lymphocytes in the bursa of Fabricius, but also many other constituent cells and tissues. They provide suitable microenvironment for differentiation and maturation of B lymphocytes. Exporting enough B cells from the bursa of Fabricius is essential for maintaining a normal peripheral B cell population and potential humoral immune responses [34]. Hence, the bursa of Fabricius as entire experimental sample was selected in this study. Indentifying genes and pathways from all cells in the bursa of Fabricius that are involved in immunomodulation in stress model can provide targets for understanding the molecular mechanisms of stress affecting immune function. Considering the influence of the heterogeneity in cell composition of the bursa of Fabricius to the interpretation of the data, the expression and function of candidate genes would be comprehensively analyzed through the annotations of database and immunohistochemical location in our subsequent studies.

Cathelicidin-B1-like (CATHB1, also known as CATHL1) was the most highly expressed gene in the bursa of Fabricius of the C_B group and the fourth highest expressed gene of B_B group. CATHB1 is one of the four cathelicidin genes encoded by the chicken genome to form a family of vertebrate-specific immune molecules produced mainly by the bursa of Fabricius[35]and showed preferential expression in bursa of Fabricius [36]. Cathelicidin genes encode host defense proteins that have both antimicrobial activity and immunomodulatory functions [37]. Diverse functions that found in mammalian cathelicidins including chemotaxis of inflammatory and immune cells, stimulation of phagocytosis and activation and differentiation of immune cells [38]. Since it is unclear in which cell the chicken cathelicidin genes are expressed, the induction of local gene expression or an influx of cathelicidin expressing cells will increase the expression of cathelicidin genes [39], the highly expression of CATHB1 need further researched.

Inflammation is one of the major reactions of the immune system against infection and irritation. Chronic inflammation is related to the expression of chemokines, cytokines and adhesion molecules. Pro-inflammatory cytokines, such as TNF-α, IL-1, IL-8, IFN-γ (interferon gamma) and IL-11, play an important role in mediating chronic inflammation. In our research, most of the cytokines in the CORT-treated group were significantly downregulated DEGs, such as IL-1 β (log2FC = − 1.426), IL-18 (log2FC = − 1.3586), IL-7 (log2FC = − 1.3558) and so on. Glucocorticoid has been shown to downregulate a variety of pro-inflammatory cytokines, including IL-1 β, IL-2, IL-4, IL-6, IL-8, IL-11, IL-12, GM-CSF and TNF-α, whereas glucocorticoid upregulated the expression of some anti-inflammatory genes including IL-1 R2(decoy receptor), CC10, IL-1 receptor antagonist and Lipocortin-1 [40]. Of these downregulated genes, the IL-6 gene was the only gene that contains a GRE promoter and recruit GR when glucocorticoids were present [41]. In the tethering mechanism, GR and transcription factors such as NF-kB and AP-1 were directly interact; thus GR can be recruited to genomic regulatory regions and can increase or decrease the transcription of genes that contain no GRE site via transactivation or trans-repression, respectively [42]. Therefore, the tethering mechanism may play an important role in the downregulation of genes that contain no GRE sites, such as IL-1 β and IL-18, when an organism is treated with CORT.

According to the results of GO term, KEGG pathway enrichment and PPI network analysis of the 1434 SDEGs, VEGFA and ILs (such as IL15) via the cytokine-cytokine receptor interaction pathway, MYD88 and TLR4 via the Toll-like receptor signaling pathway may play important roles in the regulation of immune function caused by stress. VEGF is among the known tumor-associated factors that inhibits immune cell function. It inhibits T-cell development and may contribute to tumor-induced immune suppression [43]. ILs play an important role in mediating inflammation. In this study, there was no direct interaction relationship between VEGF and IL15, which severally play an immunoregulatory role through their upstream and downstream genes. Toll-like receptor play pivotal roles in the innate immune defense mechanism [44]. MyD88 is the most common adaptor protein in the Toll-like receptor pathways. All TLRs use MyD88 alone or in combination with other adapters, except TLR3 [45]. The expression of TLR4-MyD88 in B1a cells is critical for the IgM-dependent atheroprotection that reduces the CD4+ and CD8+ T cells infiltrates and enhances the TGF-β1 expression while also reducing the lesion inflammatory cytokines TNF-α, IL-1 β, and IL-18 [46]. Astragalus polysaccharides (APS) can regulate host immune response by activating the TLR4/MYD88 signaling pathway [47]. The seizure-induced TLR4-mediated MyD88-dependent signaling pathway plays a key role in activating microglia and triggering neuron apoptosis [48]. MyD88 and TLR4 genes were significantly downregulated in our research which may indicate that it is regulated by CORT to change the immune function of the host via the TLR4/MYD88 signaling pathway.

Further research for the molecular mechanisms of stress influence immune function is still necessary to reach the goal of increased disease resistance in stressed chicken. Determining the functional impacts of MYD88 knockdown or TLR4 overexpression on B cells will verity the importance of Toll-like receptor signaling pathway. Characterization of the effects of CORT on B-cells survival and proliferation and expression of VEGFA gene would be needed to further verify our experimental result and whether they have a positive or negative effect on immune capabilities. To characterize the interaction of stress and immune function, we need further investigation the expression and functional changes in cell proliferation and migration.

Conclusion

In summary, we constructed and evaluated the stress model of chicken with administration of adding 30 mg CORT per kilogram of the basal diets. A total of 1434 SDEGs were identified in bursa of Fabricius of chicken using RNA-seq technology. Through GO enrichment, KEGG pathway and PPI network analysis, VEGFA, ILs (such as IL15) via the cytokine-cytokine receptor interaction pathway, MYD88 and TLR4 via the Toll-like receptor signaling pathway may play important roles in the regulation of immune function under stress condition with CORT administration. The results provide a basis for uncovering the molecular regulation mechanism of stress-influenced immune function in poultry.

Materials and methods

Experimental animals and tissue samples

Sixty healthy Gushi cocks of similar weight and 7 days of age were obtained from the Animal Center of Henan Agricultural University, randomly assigned to two groups, and held in a clean and relatively isolated room. All cock feeding and management were conducted in accordance with the chick feeding and management manual but without routine immunization and beak trimming in the feeding process. Control group (B_B) was fed the basal diet only, while the experimental group (C_B) was fed a basal diet, with the addition of 30 mg CORT/kg of basal diet [49](SIGMA -ALDRICH, SHANGHAI, CHINA, Purity ≥92%). Three cocks of each group were randomly selected and euthanized by exsanguination after 7 days of the dietary treatment. The stress model was evaluated by detecting the indicators related to immune/stress. The statistical analyses were performed with Graphpad Prism 5 (Graphpad Sofware, San Diego, CA), P ≤ 0.05 were considered statistically significant. The data are presented as the means ± SEM. Bursa of Fabricius tissues were harvested quickly and stored at − 80 °C until use.

RNA extraction

The total RNA was isolated from the bursa of Fabricius tissues of chickens using an RNAiso Plus kit (Takara, Kyoto, Japan). Total RNA was purified using a TruSeq RNASample Prep Kit v2 (New England Biolabs, Ipswich, MA, USA). RNA degradation and contamination were detected using 1% agarose gels. The purity of the total RNA was assessed with a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The concentration of the total RNA was measured using a Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). The integrity of the total RNA was estimated using an RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction and transcriptome sequencing

Six mRNA sequence libraries were constructed: one for every sample (B_B_1, B_B_2, B_B_3 and C_B_1, C_B_2, C_B_3). A total of 3 μg RNA per sample was needed to generate the sequencing libraries using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA), following the manufacturer’s protocol and adding index codes to attribute sequences to each sample. Briefly, mRNA with a polyA tail was purified from total RNA using poly-T oligo-attached magnetic beads by the combination of A-T. The mRNA was fragmented and the cDNA was synthesized and purified with AMPure XP system (Beckman Coulter, Beverly, USA) to allow preferential selection of cDNA fragments of 250~300 bp in length. The cDNA library was enriched by PCR with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. At last, the PCR products were purified (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system and Q-PCR. Index-coded cDNA fragment libraries were analyzed using a paired-end 2 × 150 nt Illumina HiSeq platform controlled by data collection software. The image data were outputted and transformed into raw data and stored in FASTQ (fq) format after sequencing. The Q20, Q30, and GC content of the clean data were calculated at the same time. Raw sequencing data were processed through Perl scripts to remove reads containing the adapter, reads containing more than 10% poly-N, and reads of low quality (reads containing more than 50% bases with Qphred ≤ 20). All downstream analyses were based on the resulting clean data.

Mapping and transcript abundance estimation

The chicken genome assembly (Gallus Gallus 4.0) and gene model annotation files were downloaded from Ensemble. Bowtie (version: 2.2.3) [50] was used to build the index of the reference genome, and TopHat (version: 2.1.0) [51] was used to perform genome mapping on the clean data. TopHat can generate a database of splice junctions based on the gene model annotation file, thus allowing the alignment of reads that contain gaps compared with the genome and consummating catalogs of alternative splicing events. Based on the alignment results, the sequencing data of each library were assembled separately by Cufflinks (version 2.1.1) [52]. Cufflinks were used to assemble the gene annotation files, compare them with the original gene annotation files, determine the genes that are not included in the original annotations and optimize the gene’s position, thereby supplementing and modifying the original annotation file. The transcripts were mapped to the reference annotation files using the Cuffcompare, a program in the Cufflinks package, to sort out novel genes; discover new exon regions of known genes; and optimize the starting and ending positions of known genes [53]. Then, HTSeq (version: 0.6.1) [54] was used to count the number of fragments for each gene. The FPKM (fragments per kilobase of exon model per million mapped reads) considers the effect of the sequencing depth and gene length to standardize the read count for the comparison between the different genes and samples through the use of a Perl script (FPKM = total exon fragments/mapped reads in millions × exon length in kb) [55].

Identification and functional annotation of DEGs

Differentially expressed genes of two groups (three biological replicates per group) were performed using the DESeq2 R package (version:1.14.1) by mean normalised read-counts [56]. DESeq2 provides statistical routines for determining DEGs using a model based on the negative binomial distribution. The input count values must be raw counts of sequencing reads for DESeq’s statistical model to hold, as only the actual counts allow assessing the measurement precision correctly [56, 57]. So that we did no filter of FPKM value before the DESeq analyses. The Benjamini and Hochberg’s approach were used to control the FDR of the P-value. We identified DEGs as genes with an adjusted P-value (padj) < 0.05 and FPKM > 1, and SDEGs as genes with padj < 0.05, |log2FC| ≥ 1, and FPKM > 1. Cluster analysis was used to cluster genes with the same or similar expression patterns, which might have similar functions or participate in the same biological process. Cluster analysis of a heat-map for SDEGs was performed by the pheatmap R package.

The GOseq R package [58] was used to identify significantly enriched GO terms with corrected P-values < 0.05 by SDEGs. The KEGG (http://www.genome.jp/kegg/) pathway is the major public pathway-related database for understanding the biological functions of genes, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies. We used KOBAS software to test the statistical enrichment of SDEGs in KEGG pathways [59] to identify enriched pathways and clarify the group differences in cellular pathways. Pathways with corrected P-value (q-value) < 0.05 were identified as significantly enriched by the SDEGs.

The protein-protein interaction analysis of DEGs

The PPI network of the SDEGs was analyzed using the interaction relation in the STRING protein interaction database (http://string-db.org, Organism: Chicken) by extracting the SDEGs list from the database pertaining to chickens and selected the interaction relation with confidence score > 0.7. The PPI networks were visualized by Cytoscape (version 3.6.1). Based on the scale-free property of interaction networks, hub proteins were found by counting the number of interactions of each network node has with other nodes [60]. Network module analysis was performed through MCODE in Cytoscape with default parameters (Degree Cutoff: 2, Node Score Cutoff: 0.2, K-Core: 2, Max. Depth: 100).

qRT-PCR

We randomly selected 12 genes (including 4 upregulated genes and 8 downregulated genes) for qRT-PCR in the two groups to confirm the reproducibility and accuracy of the RNA-Seq gene expression data. Using an RNAiso Plus kit (Takara, Kyoto, Japan), total RNA was isolated from bursa of Fabricius tissues of chickens. After the check of RNA quality, as described in the section on RNA extraction, reverse transcription was performed using the PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara, Kyoto, Japan) according to the manufacturer’s protocol.

The qRT-PCR experiments were performed using a LightCycler® 96 Real-Time PCR system (Roche Applied Science) in a 25-μL reaction volume containing 12.5 μL of 2 × SYBR® Premix Ex TaqTM II (Tli RNaseH Plus) (Takara, Kyoto, Japan), 1.25 μL each of the forward and reverse primers (10 μM), 8 μL of deionized water, and 2 μL (approximately 100 ng) of cDNA. The GAPDH gene was used as the reference gene, and the primers of 12 genes were designed using the Primer-BLAST website of NCBI (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primers are shown in the Additional file 14: Table S10. The thermal cycling conditions were 3 min at 95 °C, followed by 37 reaction cycles (95 °C for 30 s, 60 °C for 30 s, and 72 °C for 30 s), and an extension for 10 min at 72 °C. We calculated the relative gene expression levels with the comparative CT method (also referred to as the 2-△△CT method) [61], with three replicates each reaction. The data are presented as the means ± SEM.

Abbreviations

ACTH:

Adrenocorticotropic hormone

AP-1:

Activator protein-1

B2M:

β-2-microglobulin

CATHB1:

Cathelicidin-B1-like

CORT:

Corticosterone

CRH:

Corticotropin-releasing hormone

GR:

Glucocorticoid receptor

GRE:

Glucocorticoid response element

IGLL1:

Immunoglobulin lambda-like polypeptide 1

MR:

Mineralocorticoid receptor

MYD88:

Myeloid differentiation primary-response gene 88

NF-Κb:

Nuclear factor-kappa B

RSAD2:

Radical S-adenosyl methionine domain-containing 2

VEGFA:

Vascular endothelial growth factor

References

  1. Chrousos GP. Stress and disorders of the stress system. Nat Rev Endocrinol. 2009;5(7):374.

    Article  CAS  PubMed  Google Scholar 

  2. Kumar B, Manuja A, Aich P. Stress and its impact on farm animals. Frontiers in Bioscience. 2012;4(9):1759.

    Article  Google Scholar 

  3. Sun GR, Yan LI, Kang XT, Tian YD, Zhang H, Kui LI. Effect of beak trimming on the apoptosis and related protein expression of chicken Thymus. Chinese Journal of Animal & Veterinary Sciences. 2011;42(7):1000–1006.

  4. Kaiser P, Wu Z, Rothwell L, Fife M, Gibson M, Poh TY, Shini A, Bryden W, Shini S. Prospects for understanding immune-endocrine interactions in the chicken. Gen Comp Endocrinol. 2009;163(1–2):83–91.

    Article  CAS  PubMed  Google Scholar 

  5. Shini S, Huff GR, Shini A, Kaiser P. Understanding stress-induced immunosuppression: exploration of cytokine and chemokine gene profiles in chicken peripheral leukocytes. Poult Sci. 2010;115(1):117–23.

    Google Scholar 

  6. Wu CC, Dorairajan T, Lin TL. Effect of ascorbic acid supplementation on the immune response of chickens vaccinated and challenged with infectious bursal disease virus. Vet Immunol Immunopathol. 2000;74(1–2):145–52.

    Article  CAS  PubMed  Google Scholar 

  7. Sarkar P, Ghosh S, Batabyal S, Chatterjee S. Biochemical stress responses of broiler chickens during transport. Indian Journal of Animal Research. 2013;47(1):29–34.

    Google Scholar 

  8. Fisinin VI, Miftakhutdinov AV, Anosov DE. Pharmacological prevention of stress during chicken debeaking. Russ Agric Sci. 2016;42(1):97–100.

    Article  Google Scholar 

  9. Khatoon ZAA. Heat stress in poultry and the beneficial effects of ascorbic acid (vitamin C) supplementation during periods of heat stress. Worlds Poultry Science Journal. 2013;69(1):135–52.

    Article  Google Scholar 

  10. Tang J, Chen Z. The protective effect of γ-aminobutyric acid on the development of immune function in chickens under heat stress. Journal of animal physiology. Animal Nutrition. 2016;100(4):768–77.

    CAS  Google Scholar 

  11. Zhao JP, Lin H, Jiao HC, Song ZG. Corticosterone suppresses insulin- and NO-stimulated muscle glucose uptake in broiler chickens ( Gallus Gallus domesticus ). Comp Biochem Physiol C Toxicol Pharmacol. 2009;149(3):448–54.

    Article  CAS  PubMed  Google Scholar 

  12. Duan Y, Fu W, Wang S, Ni Y, Zhao R. Cholesterol deregulation induced by chronic corticosterone (CORT) stress in pectoralis major of broiler chickens. Comparative Biochemistry & Physiology Part a Molecular & integrative. Physiology. 2014;176:59–64.

    CAS  Google Scholar 

  13. Dong H, Lin H, Jiao HC, Song ZG, Zhao JP, Jiang KJ. Altered development and protein metabolism in skeletal muscles of broiler chickens (Gallus gallus domesticus) by corticosterone. Comp Biochem Physiol A Mol Integr Physiol. 2007;147(1):189–95.

    Article  CAS  PubMed  Google Scholar 

  14. Gao J, Lin H, Song ZG, Jiao HC. Corticosterone alters meat quality by changing pre-and Postslaughter muscle metabolism. Poult Sci. 2008;87(8):1609.

    Article  CAS  PubMed  Google Scholar 

  15. Berger S, Romero LM, Kalko EK, Vitousek MN, Rodl T. Corticosterone suppresses immune activity in territorial Galapagos marine iguanas during reproduction. Hormones & Behavior. 2005;47(4):419–29.

    Article  CAS  Google Scholar 

  16. Lin H, Decuypere E, Buyse J. Oxidative stress induced by corticosterone administration in broiler chickens (Gallus gallus domesticus) 1. Chronic exposure. Comparative Biochemistry & Physiology Part B Biochemistry. Mol Biol. 2004;139(4):737–44.

    CAS  Google Scholar 

  17. Cui YL, Zheng SX, Hu YF, Peng YX. Effect of heat-stress on the immune organ of chickens. In: Journal of Agricultural University of Hebei; 2004.

    Google Scholar 

  18. Wang L, Luo ZG, Cai MC, Fu PH, Zhou P, Zuo FY. Effect of heat stress on serum ionic Concentration,Enzymatic Activity,Antioxidant and immune index in beef cattle. In: Journal of Southwest China Normal University; 2017.

    Google Scholar 

  19. Ju XH, Yong YH, He JC, Deng FC: Effect of Heat Stress on Baseline Immune and Blood Biochemical Index in Ba-Ma Miniature Pig. Chinese Journal of Animal Science. 2009;45(13):51–54.

  20. Lin W, Li L, Chen J, Li D, Hou J, Guo H, Shen J. Long-term crowding stress causes compromised nonspecific immunity and increases apoptosis of spleen in grass carp (Ctenopharyngodon idella). Fish & Shellfish Immunology. 2018;(28):540–45.

  21. Xiang Y, Yan H, Zhou J, Zhang Q, Hanley G, Caudle Y, Lesage G, Zhang X, Yin D. The role of toll-like receptor 9 in chronic stress-induced apoptosis in macrophage. PLoS One. 2015;10(4):e0123447.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Zhang Y, Li D, Han R, Wang Y, Li G, Liu X, Tian Y, Kang X, Li Z. Transcriptome analysis of the pectoral muscles of local chickens and commercial broilers using Ribo-zero ribonucleic acid sequencing. PLoS One. 2017;12(9):e0184115.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Li Z, Liu X, Zhang P, Han R, Sun G, Jiang R, Wang Y, Liu X, Li W, Kang X. Comparative transcriptome analysis of hypothalamus-regulated feed intake induced by exogenous visfatin in chicks. BMC Genomics. 2018;19(1):249.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Zhang M, Li DH, Li F, Sun JW, Jiang RR, Li ZJ, Han RL, Li GX, Liu XJ, Kang XT. Integrated analysis of MiRNA and genes associated with meat quality reveals that Gga-MiR-140-5p affects intramuscular fat deposition in chickens. Cellular Physiology & Biochemistry International Journal of experimental cellular physiology biochemistry & Pharmacology. 2018;46(6):2421.

    Article  CAS  Google Scholar 

  25. Bao Y, Liang C, Sa R, Zhong R, Xing H, Zhang H. Transcriptome profile analysis of breast muscle tissues from high or low levels of atmospheric Ammonia exposed broilers (Gallus gallus). PLoS One. 2016;11(9):e0162631.

    Article  CAS  Google Scholar 

  26. Sun H, Jiang R, Xu S, Zhang Z, Xu G, Zheng J, Qu L. Transcriptome responses to heat stress in hypothalamus of a meat-type chicken. Journal of Animal Science and Biotechnology. 2015;6(3):6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Paramithiotis E, Michael JHR. B cell emigration directly from the cortex of lymphoid follicles in the bursa of Fabricius. Eur J Immunol. 1994;24(2):458–63.

    Article  CAS  PubMed  Google Scholar 

  28. Babacanoğlu E, Yalçin S, Uysal S. Evaluation of a stress model induced by dietary corticosterone supplementation in broiler breeders: effects on egg yolk corticosterone concentration and biochemical blood parameters. Br Poult Sci. 2013;54(6):677–85.

    Article  PubMed  CAS  Google Scholar 

  29. Dostert A, Heinzel T. Negative glucocorticoid receptor response elements and their role in glucocorticoid action. Curr Pharm Des. 2004;10(23):2807–16.

    Article  CAS  PubMed  Google Scholar 

  30. He Y, Xu Y, Zhang C, Gao X, Dykema KJ, Martin KR, Ke J, Hudson EA, Khoo SK, Resau JH. Identification of a lysosomal pathway that modulates glucocorticoid signaling and the inflammatory response. Sci Signal. 2011;4(180):ra44.

  31. Luecke HF, Yamamoto KR. The glucocorticoid receptor blocks P-TEFb recruitment by NFkappaB to effect promoter-specific transcriptional repression. Genes Dev. 2005;19(9):1116.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Vahdatpour T, Adl KN, Nezhad YE, Sis NM, Riyazi SR, Vahdatpour S. Effects of corticosterone intake as stress-alternative hormone on broiler chickens: performance and blood parameters. Asian Journal of Animal & Veterinary Advances. 2009;4(1):16–21.

    Article  CAS  Google Scholar 

  33. Xie WY, Hou XY, Yan FB, Sun GR, Han RL, Kang XT. Effect of γ-aminobutyric acid on growth performance and immune function in chicks under beak trimming stress. Anim Sci J. 2013;84(2):121–9.

    Article  CAS  PubMed  Google Scholar 

  34. Monson MS, Goor AGV, Ashwell CM, Persia ME, Rothschild MF, Schmidt CJ, Lamont SJ. Immunomodulatory effects of heat stress and lipopolysaccharide on the bursal transcriptome in two distinct chicken lines. BMC Genomics. 2018;19(1):643.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. Sang IL, Jang HJ, Jeon M, Mi OL, Kim JS, Jeon IS, Byun SJ. Transcriptional regulation of cathelicidin genes in chicken bone marrow cells. Poult Sci. 2016;95(4):912.

    Article  CAS  Google Scholar 

  36. Lee MO, Jang HJ, Rengaraj D, Yang SY, Han JY, Lamont SJ, Womack JE. Tissue expression and antibacterial activity of host defense peptides in chicken. BMC Vet Res. 2016;12(1):231.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Goitsuka R, Chen CL, Benyon L, Asano Y, Kitamura D, Cooper MD. Chicken cathelicidin-B1, an antimicrobial guardian at the mucosal M cell gateway. Proc Natl Acad Sci U S A. 2007;104(38):15063.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zanetti M. The role of cathelicidins in the innate host defenses of mammals. Current Issues in Molecular Biology. 2005;7(2):179–96.

    CAS  PubMed  Google Scholar 

  39. Rodriguez-Lecompte JC, Yitbarek A, Cuperus T, Echeverry H, Van DA. The immunomodulatory effect of vitamin D in chickens is dose-dependent and influenced by calcium and phosphorus levels. Poult Sci. 2016;95(11):2547–56.

    Article  CAS  PubMed  Google Scholar 

  40. Barnes PJ. Anti-inflammatory actions of glucocorticoids: molecular mechanisms. Clin Sci. 1998;94(6):557–72.

    Article  CAS  Google Scholar 

  41. Ray A, Laforge KS, Sehgal PB. On the mechanism for efficient repression of the interleukin-6 promoter by glucocorticoids: enhancer, TATA box, and RNA start site (Inr motif) occlusion. Molecular & Cellular Biology. 1990;10(11):5736–46.

    Article  CAS  Google Scholar 

  42. Muzikar KA, Nickols NG, Dervan PB. Repression of DNA-binding dependent glucocorticoid receptor-mediated gene expression. Proc Natl Acad Sci U S A. 2009;106(39):16598–603.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Ohm JE, Gabrilovich DI, Sempowski GD, Kisseleva E, Parman KS, Nadaf S, Carbone DP. VEGF inhibits T-cell development and may contribute to tumor-induced immune suppression. Blood. 2003;101(12):4878–86.

    Article  CAS  PubMed  Google Scholar 

  44. Wang M, Wang L, Jia Z, Yi Q, Song L. The various components implied the diversified toll-like receptor (TLR) signaling pathway in mollusk Chlamys farreri. Fish Shellfish Immunol. 2018;74:205–12.

    Article  CAS  PubMed  Google Scholar 

  45. Naro C, Sette C. Dissecting a hub for immune response: modeling the structure of MyD88. Structure. 2016;24(3):349–51.

    Article  CAS  PubMed  Google Scholar 

  46. Hosseini H, Li Y, Kanellakis P, Tay C, Cao A, Liu E, Peter K, Tipping P, Toh BH, Bobik A. Toll-like receptor (TLR)4 and MyD88 are essential for Atheroprotection by peritoneal B1a B cells. J Am Heart Assoc. 2016;5(11):e002947.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Zhou L, Liu Z, Wang Z, Yu S, Long T, Zhou X, Bao Y. Astragalus polysaccharides exerts immunomodulatory effects via TLR4-mediated MyD88-dependent signaling pathway in vitro and in vivo. Sci Rep. 2017;7:44822.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Hu QP, Mao DA. Histone deacetylase inhibitor SAHA attenuates post-seizure hippocampal microglia TLR4/MYD88 signaling and inhibits TLR4 gene expression via histone acetylation. BMC Neurosci. 2016;17(1):22.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Malheiros R, Moraes V, Collin A, Decuypere E, Buyse J. Free diet selection by broilers as influenced by dietary macronutrient ratio and corticosterone supplementation. 1. Diet selection, organ weights, and plasma metabolites. Poult Sci. 2003;82(1):123–31.

    Article  CAS  PubMed  Google Scholar 

  50. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7(3):562.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.

    Article  CAS  PubMed  Google Scholar 

  55. Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621.

    Article  CAS  PubMed  Google Scholar 

  56. Anders S, Huber W: Differential expression of RNA-Seq data at the gene level–the DESeq package. Heidelberg, Germany: European Molecular Biology Laboratory (EMBL) 2012.

  57. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(2):R14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    Article  CAS  PubMed  Google Scholar 

  60. Xu H, Ma J, Wu J, Chen L, Sun F, Qu C, Zheng D, Xu S. Gene expression profiling analysis of lung adenocarcinoma. Brazilian Journal of Medical & Biological Research. 2016;49(3):S0100.

    Article  CAS  Google Scholar 

  61. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nature Protocol. 2008;3(6):1101.

    Article  CAS  Google Scholar 

Download references

Acknowledgments

We are grateful to the editor and reviewers for their comments and suggestions.

Funding

This research was supported by the Program for Innovation Research Team of Ministry of Education (IRT16R23), Natural Science Foundation of Henan Province (162300410141), Key Scientific Research Project Plan of Henan Province Colleges and Universities (14A230018, 17A230015), Science and Technology Planning Project of Henan Province (162300410172).

Availability of data and materials

All raw sequence data generated in this study were deposited in the NCBI Sequence Read Archive (https://submit.ncbi.nlm.nih.gov/subs/sra/) under accession number SRR7959183–7959188 (BioProject number PRJNA494531; https://www.ncbi.nlm.nih.gov/bioproject/494531). All other supporting data are included as additional files.

Author information

Authors and Affiliations

  1. College of Animal Science and Veterinary Medicine, Henan Agricultural University, Zhengzhou, 450002, China

    Yanhua Zhang, Yanting Zhou, Guirong Sun, Kui Li, Zhuanjian Li, Aru Su, Xiaojun Liu, Guoxi Li, Ruirui Jiang, Ruili Han, Yadong Tian, Xiangtao Kang & Fengbin Yan

Authors
  1. Yanhua Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yanting Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Guirong Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Kui Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zhuanjian Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Aru Su
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Xiaojun Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Guoxi Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Ruirui Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Ruili Han
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Yadong Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Xiangtao Kang
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Fengbin Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YZha performed the data analysis and wrote the main manuscript; FY conceived the study, designed the experiments, and helped to draft the manuscript; YZho performed the RNA-related experiments; XK and XL participated in the experiment design and the discussion; RH, GS and KL participated in the GO analysis and prepared the GO-related figures and tables; YT and GL performed the pathway analysis and prepared the pathway-related figures and tables; AS and RJ contributed to the bioinformatics analysis of the sequencing data; ZL participated in the gene function analysis. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Fengbin Yan.

Ethics declarations

Ethics approval and consent to participate

The experiments and animal care were performed in accordance with the protocol approved by the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, 2004), and the protocol was approved by the Institutional Animal Care and Use Committee (Use the Animal Care Committee of the College of Animal Science and Veterinary Medicine, Henan Agricultural University, China; Permit Number: 17–0118). All experiments and methods were performed according to approved guidelines to minimize animal suffering.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Characteristics of the reads from bursa of Fabricius libraries obtained from 2 groups. 1Multiple mapped = number of clean reads and the ratio that matched two or more positions in the genome. 2Uniquely mapped = number of clean reads and the ratio that matched only one position in the genome. (DOCX 17 kb)

Additional file 2:

Figure S1. Pearson correlation between samples. (TIF 362 kb)

Additional file 3:

Figure S2. Venn diagram of global genes expressed in two groups. Red indicates the genes specific expressed in experimental group, blue indicates the genes specific expressed in the control group and purple indicates the genes expressed in both group. (TIF 504 kb)

Additional file 4:

Table S2. The gene expression level of the two groups. Sheet 1 showed the top 10 highest expression genes in two groups. Sheet 2 showed the FPKM and read count of all genes. Sheet 3 showed the mean normalised read-counts by DESeq2. (XLSX 3969 kb)

Additional file 5:

Table S3. The DEGs in the C_B group compared with the B_B group. Sheet 1 showed the SDEGs, sheet 2 showed the downregulated SDEGs (padj < 0.05, |log2FC| ≥ 1, and FPKM > 1), sheet 3 showed the upregulated SDEGs, sheet 4 showed the DEGs (padj < 0.05, FPKM > 1). (XLSX 1358 kb)

Additional file 6:

Table S4. The GO enrichment of SDEGs of the C_B group and B_B group. Sheet 1 showed the top 30 most significantly enriched GO terms and sheet 2 showed all the enriched GO terms. (XLSX 974 kb)

Additional file 7:

Table S5. The KEGG enrichments of the differentially expressed genes of the C_B group and B_B group. Sheet 1 showed the top 20 most significantly enriched KEGG pathways and sheet 2 showed all the enriched KEGG pathways. (XLSX 47 kb)

Additional file 8:

Table S6. The correlation coefficient between receptor genes expression and ligand genes expression in the cytokine-cytokine receptor. (XLSX 21 kb)

Additional file 9:

Table S7. Protein-protein interaction network for all SDEGs. Sheet 1 showed the interactions with a confidence score > 0.7, sheet 2 showed the log2FoldChange values of genes in sheet1, sheet 3 showed the KEGG analysis result of the genes in the 19 modules. (XLSX 29 kb)

Additional file 10:

Figure S3. Protein-protein interaction network for all SDEGs. The color of nodes indicates the degree of foldchange, the value was “-log2FoldChange”. Red color of label indicates proteins which exceeded 8 interactions with others, while black color indicates proteins which below 8 interactions. The edge between proteins represent interactions between them and the confidence score were more than 0.7. (TIF 695 kb)

Additional file 11:

Figure S4. 19 modules of the protein-protein interaction network with node > 3 and a MCODE score > 3. The color of nodes indicates the degree of foldchange, the value was “-log2FoldChange”. The edge between proteins represent interactions between them and the confidence score were more than 0.7. (TIF 617 kb)

Additional file 12:

Table S8. Effects of CORT treatment on the serum biochemical indexes of chickens. Data are shown as the mean ± SE. Different lowercase letters (a and b) in same column indicate significant differences among the C_B and B_B groups (P < 0.05). (DOCX 13 kb)

Additional file 13:

Table S9. Effects of CORT treatment on immune related indexes of chickens. Data are shown as the mean ± SE. Different lowercase letters (a and b) in same column indicate significant differences among the C_B and B_B groups (P < 0.05). (DOCX 15 kb)

Additional file 14:

Table S10. List of the genes and primers used for qRT-PCR validation. (XLSX 12 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, Y., Zhou, Y., Sun, G. et al. Transcriptome profile in bursa of Fabricius reveals potential mode for stress-influenced immune function in chicken stress model. BMC Genomics 19, 918 (2018). https://doi.org/10.1186/s12864-018-5333-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5333-2

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us