Transcriptome analysis revealed the roles of long non-coding RNA and mRNA in the bursa of Fabricius during pigeon (Columba livia) development

The bursa of Fabricius (BF) is the critical humoral immune organ to birds, playing an essential role in B lymphocyte differentiation. However, unlike other poultries, surgical removal of pigeon BF did not limit humoral immune responsiveness. To investigate the expression profiles and the potential role of mRNA and long non-coding RNA (LncRNA) in squab BFs, transcriptome analysis was performed by RNA-Sequencing (RNA-Seq) over three developmental stages (1-day, 13 and 26 days old). We identified 13,072 mRNAs and 19,129 lncRNAs, of which 2,752 mRNAs and 1,515 lncRNAs were differential expressed (DE) in pigeon BFs over three developmental stages. Cluster analysis presented different expression patterns in DE mRNAs and lncRNAs. Functional enrichment analysis revealed that DE lncRNAs and mRNAs with distinct expression patterns might play crucial roles in the immune system process and tissue morphogenesis. In particular, some DE genes and lncRNAs with higher expression levels in 13D or 26D are related to lymphocyte activation and differentiation, adaptive immune response, positive regulation of immune response, leukocyte migration, etc. Protein-protein interaction (PPI) network and Molecular Complex Detection (MCODE) analysis sreened six significant modules containing 37 genes from immune-related DE gene cluster, which is closely linked in B cell activation, lymphocyte differentiation, B cell receptor signaling pathway, etc. Our study characterizes mRNA and lncRNA transcriptomic variability in pigeon BFs over different developmental stages and enhances understanding of the mechanisms underlying physiological functions of pigeon BF.

Introduction

Mammals and birds evolved from a common reptilian ancestor and have many common immunological systems (1). Meanwhile, they also developed some distinctly different characteristics. For instance, unlike human beings and mice, bird B-cells develop in the bursa of Fabricius (2). The bursa of Fabricius, an obscure sac-like structure attached to the proctodeal region of the cloaca (1), is the central humoral immune organ unique to birds (3). BF is responsible for the amplification and differentiation of B lymphoid progenitors within its follicular microenvironment (4). Fully differentiated B lymphocytes migrate from BF to peripheral lymphoid organs to colonize, reproduce and perform important immune functions (5). Of note, BF contains several biologically-active factors, including bursin, brusal peptide 11 (BP11), bursal pentapeptide I (BPP-I), etc. These factors are involved in B cell differentiation, regulation of B-cell development, antigen-specific immune responses (6), antiproliferation of tumor cell, and immunomodulatory activity (7). Due to the importance of BF in B-cell development and immunoglobulin production, surgical removal or pathological changes in BF would result in reduction of humoral immunocompetence (8), and an increase in the susceptibility to infections at certain birds (9). Additionally, BF size and weight might influence the immunity of avian (10, 11).

In China, the pigeon (Columba livia) is an indispensable agricultural economic animal for meat and eggs, and has become the fourth-largest domestic poultry by breeding scale, following chicken, duck and goose. Pigeon squabs are altrices, and have an extraordinarily high growth rate. Pigeon BF during development undergoes hypertrophia and hyperplasia of lobules constituting the plicae, and the gradual and continuous increasing of lymphocytes (12). After treatment with cyclophosphamide, pigeon BF produced a rapid atrophy and an intense depletion of lymphoid cells (13). However, unlike other poultries (chickens, White Pekin ducks, Japanese quail, etc.), bursectomy of pigeon in early neonatal period did not affect the development of humoral immune responsiveness (14). This phenomenon indicates that the pigeon BF development process may be different from other poultries in some way.

Chicken BF transcriptome research across distinct developmental stages indicated that DE genes are related in defence response and some signal pathways (Jak-STAT, TLR, Wnt and MAPK, etc) might be involving in the regulation of BF development process (2). Additionally, long non-coding RNAs (lncRNAs) have been found to play crucial roles during organ development in recent years (15). lncRNAs are defined as non-protein coding transcripts longer than 200 nucleotides (16), which regulate gene expression at the epigenetic, transcriptional and post-transcriptional levels via different mechanisms (17). To date, lncRNAs and mRNAs closely linked to the functional development of pigeon BF have not been identified. Here, we aim to explore the lncRNA and mRNA expression profiles in BF and their potential roles during pigeon development by RNA sequencing technology. Our findings will provide insights into the transcriptional variations in pigeon BF development and enhance the understanding physiological functions of squab BF.

Results

Phenotypic measurements

In this study, we investigated morphological changes, BF weight, and BF index across three different development stages of pigeons. As depicted in Figure 1A– C, BF size and weight conspicuously increased from 1-day-old to 26-day-old. BF indexes in 13-day-old and 26-day-old squabs are conspicuously higher than those of 1-day-old pigeons. However, BF indexes did not exhibit a significant difference between 13-day-old and 26-day-old pigeon squabs. From the results of BF morphological changes, lymphoid follicles in BFs gradually enlarged from 1-day-old to 26-day-old (Figure 1D– F).

FIGURE 1

Figure 1 BF weight, BF index and histomorphological changes at the different development stages of pigeon squabs. (A).BF size (B). BF weight (C). BF index (n = 13) (D-F). Histomorphological characteristics in pigeon BF over the different development stages. ** indicates P < 0.01.

Transcriptome sequencing and assembly

To systematically investigate lncRNA and mRNA expression in BF during pigeon development, we yielded a total of 103.56 Gb raw sequence data from 9 libraries (each stage constructed 3 libraries). After quality control, 340.91 million clean reads (102.27 Gb data) were retained for further analysis. Approximately 81.51~94.02% of the clean data were mapped to the pigeon genome (Supplementary Table 2).

In total, 13,072 mRNAs and 19,129 lncRNAs (containing 2,540 known and 16,589 novel lncRNAs) were identified in pigeon BFs during development (Supplementary Table 3). Based on characteristics comparisons with mRNAs, lncRNAs have lower coding potentials, fewer exon numbers, lower expression levels and shorter transcript lengths (Supplementary Figure 1). We screened out 2,752 DE mRNAs and 1,515 DE lncRNAs (Supplementary Table 4), accounting for 21.05% and 7.92% of total identified mRNAs and lncRNAs, respectively. More specifically, three pairwise comparisons between the stages (13D vs. 1D, 26D vs. 1D, 26D vs. 13D) were performed to identify differentially expressed mRNAs and lncRNAs, and obtained 2354, 1891, 15 DE mRNAs (Figure 2A) and 878, 1084, 11 DE lncRNAs (Figure 2B), respectively. Both differentially expressed mRNAs and lncRNAs numbers between 1D and elder age stage (13D or 26D) are much greater than those between 13D and 26D. Furthermore, three DE genes (AKR1D1, VSIG1 and BCL2A1) overlapped in all three comparisons (Figure 2C).

FIGURE 2

Figure 2 Differentially expressed mRNAs and lncRNAs during squab BF development. (A–B). Volcano plot of the DE mRNAs(A) and lncRNAs (B) between two different development stages. The x-axis indicates the difference in expression level on a log2 (fold change). The y-axis represents the corresponding false discovery rate on a negative log 10(FDR). (C). Venn plot shows the number of overlapping DE mRNAs and lncRNAs in different developmental stages.

Subsequently, hierarchical clustering analyses were conducted on the expression profiles of differentially expressed mRNAs and lncRNAs. As depicted in Figure 3A–B, DE mRNAs and lncRNAs expression data in pigeon BFs clustered into two groups according to developmental stages, respectively. The elder squabs (13D and 26D) were clustered into a subgroup, and separated from 1-day-old pigeon squabs. It implied that the developmental time leads to the different mRNA and lncRNA transcriptomes in pigeon BFs.

FIGURE 3

Figure 3 Hierarchical clustering analysis of differentially expressed mRNAs and lncRNAs in pigeon BFs. Heatmap representing the Z-score for normalized expression values of DE mRNAs (A) and lncRNAs (B). The heatmap is drawn with Pretty Heatmap at ImageGP.

Expression patterns and potential functions of DE mRNAs during BF development

The expression patterns of DE mRNAs in different BF developmental stages were analyzed by using the TCseq R package. DE mRNAs were divided into three categories (Figure 4A–B, Supplementary Table 5). Cluster 1 comprised 1305 mRNAs and their expression levels in 26D are lower than those in 1D and 13D. Cluster 2 containing 337 mRNAs showed higher expression in 13D or 26D compared to 1D, while the expression of 1110 mRNAs in cluster 3 was rapidly downward from 1D to 13D and then increased or slowly decreased in 26D (Figure 4B).

FIGURE 4

Figure 4 Distinct expression patterns of differentially expressed mRNAs across the BF development stages. (A). The heatmap shows the expression patterns of the three clusters. Each row represents a mRNA. The color indicates the Z-score of normalized expression values. (B). Cluster analysis of the DE mRNA expression patterns across the BF development stages. Membership values indicate the goodness of fit for mRNAs in a particular cluster.

To explore the potential functions of the mRNAs in different clusters, functional enrichment analyses were conducted with Metascape. As depicted in Figure 5, genes in cluster 1 and cluster 3 were mainly enriched for the GO terms related to tissue morphogenesis and development, whereas those in cluster 2 were enriched in the immune related GO terms, such as lymphocyte activation (GO:0046649), B cell activation (GO:0042113), lymphocyte differentiation (GO:0030098), leukocyte differentiation (GO:0002521), etc (Supplementary Table 6). KEGG pathway enrichment analysis also showed that enriched pathways for genes in cluster 2 are also linked with immune, including Cytokine-cytokine receptor interaction (ko04060), Primary immunodeficiency (ko05340), T cell receptor signaling pathway(ko04660), B cell receptor signaling pathway(ko04662), etc (Figure 6, Supplementary Table 6).

FIGURE 5

Figure 5 Function enrichment analysis of DE mRNAs in different clusters. The x-axis indicates –log (P - value). The y-axis indicates functional categories.

FIGURE 6

Figure 6 KEGG pathway analysis for DE mRNAs in different clusters. The x-axis indicates –log (P - value). The y-axis indicates KEGG pathways.

PPI network and MCODE enrichment analysis

Given genes in cluster 2 are closely related to immune, these genes were employed to construct PPI network. (Figure 7A). The PPI network of cluster 2 consisting of 179 nodes and 315 edges was constructed in the different protein interaction databases including STRING (18), BioGrid (19), OmniPath (20), InWeb_IM (21). MCODE algorithm was employed to identify highly connected network components from the PPI network. Six significant modules containing 37 genes were screened from our constructed PPI network by MCODE (Figure 7B, Supplementary Table 7). These significant vital modules showed functions including B cell activation, lymphocyte differentiation, positive regulation of leukocyte cell-cell adhesion, DNA repair, ion homeostasis, etc (Supplementary Table 7).

FIGURE 7

Figure 7 PPI interaction networks and core modules for DE genes in cluster 2. (A). PPI interaction network. (B). Modules of PPI.

Expression patterns and potential functions of DE lncRNAs during BF development

Like DE mRNAs, the expression patterns of DE lncRNAs were also divided into three categories (Figure 8A, B, Supplementary Table 5). Cluster 1 comprising 355 lncRNAs showed higher expression level in 13D or 26D compared to 1D. Cluster 2 contains 469 lncRNAs, and their expression levels in 26D are lower than those in 1D and 13D. Expression of 691 lncRNAs in cluster 3 was rapidly downward from 1D to 13D and then increased or slowly decreased in 26D (Figure 8B). To explore the potential functions of the lncRNAs in different clusters, cis target genes of DE lncRNAs in each cluster were predicted and subjected to functional analysis. As a result, 21, 72 and 129 potential target genes were detected in cluster 1-3, respectively (Supplementary Table 8). Target genes of DE lncRNA in cluster 1 are involved in developmental growth (GO:0048589), leukocyte migration (GO:0050900) and positive regulation of locomotion (GO:0040017), etc. Those in cluster 2 and 3 are mainly enriched for GO terms associated with development and morphogenesis (Figure 8C and Supplementary Table 9).

FIGURE 8

Figure 8 Distinct expression patterns of differentially expressed lncRNAs across the BF development stages. (A). The heatmap shows the expression patterns of the 3 clusters. Each row represents a lncRNA. The color indicates the Z-score of normalized expression values. (B). Cluster analysis of the DE lncRNA expression patterns across the BF development stages. Membership values indicate the goodness of fit for lncRNAs in a particular cluster. (C). Representative enrichment categories of target genes of DE lncRNA in different clusters.

qPCR validation

To confirm the RNA-seq results, we randomly selected eight DE genes and lncRNAs to conduct qPCR assay using three independent samples. The results were in line with our sequencing result (Pearson r =0.960 ± 0.073) (Figure 9), which highlighted the reliability of our sequencing data.

FIGURE 9

Figure 9 Validation of the sequencing data using qPCR. Four randomly selected protein coding genes and four lncRNAs in pigeon BFs were validated by real-time qPCR (n = 3). The Pearson product–moment correlation coefficient (r) was calculated using R.

Discussion

BF is acknowledged as a central humoral immune organ unique to poultries. BF development goes through three phases, including the rapid growth phase (first 3 weeks after hatch), the plateau phase, and the regression phase (22). Three development stages (1d, 13d, and 26d) were selected in this study, representing newly hatched, rapid development and plateau stage, respectively. We identified 13,072 mRNAs and 19,129 lncRNAs from three development phases of squab BFs. BF development process is accompanied by changes in gene expression. Specifically, 2,752 mRNAs and 1,515 lncRNAs were identified as differentially expressed. Moreover, it was found that DE mRNAs and lncRNAs numbers between 1D and elder age are much greater than those between 13D and 26D. It might be associated with greater BF size or weight difference between 1D and 13D or 26D.

Generally, the gene expression pattern is well connected with gene function (23). In our study, gene expression pattern analysis found that identified DE protein coding genes could be clustered into three groups. GO functional enrichment analysis found genes in cluster 1 and cluster 3 were implicated in the tissue morphogenesis and development. Of special interest and importance is DE genes in another cluster (cluster 2) which have higher expression level on 13D or 26D compared to 1D. Their functional annotations are principally involved in immune system process-related GO categories, including lymphocyte activation, B cell activation, lymphocyte differentiation, leukocyte differentiation and T cell differentiation, etc. It has been well established that BF is the location of B lymphocyte differentiation and maturation in avians. BF provides an environment to support the differentiation of B lineage cells (24), consistent with our functional enrichment analysis results. Furthermore, KEGG pathway enrichment analysis revealed that these DE genes in cluster2 are involved in Cytokine-cytokine receptor interaction, B cell receptor signaling pathway and MAPK signaling pathway, which have been demonstrated to be relevant to BF development and functions (2, 25). Representative DE genes involved in immune system process contain CD79b, BCL6, CD83, CD74, etc. Among these, CD79b also known as B29 and Ig beta, is an essential component of the B cell receptor. CD79b exhibited strict B cell-specific expression (26) and expressed from early in B-cell development to the plasma cell stage (27). CD79b plays a critical role in B cell development. Targeted knockout of CD79b gene severely disrupts pre-B-cell development (28, 29). During chicken BF development, CD79b gene is also identified as differentially expressed and shows a similar expression tendency with our results (2). Bcl-6 is not only involved in the regulation of effector and memory differentiation of B lymphocytes, but also influences effector and memory differentiation in CD4+ T cells and CD8+ T cells (30). CD83 participates in the regulation of T- and B-lymphocyte maturation and in the regulation of their peripheral responses (31). CD74 also exerts modulation functions in the process of T-cell and B-cell developments (32). These critical genes in BF exhibit higher expression levels at 13D and 26D compared with 1D, implying BF with larger size might play more positive roles in lymphocyte differentiation and development. Several lines of evidence also support this view (11, 33).

lncRNAs are crucial regulators of cell differentiation and development (34), which participate in the regulation of gene expression at multiple levels. By interacting with RNA, chromatin, and protein, lncRNAs modulate mRNA stability, chromatin structure, and the function of proteins (35). lncRNAs have the potential to regulate immune cell function and development. For instance, Tiago F. Brazão and colleagues identified 4516 lncRNAs expressed in 11 stages of mice B-cell development and activation (36). Our results found that in addition to some protein-coding genes related immune system process and tissue morphogenesis, 1,515 lncRNAs are identified as differentially expressed during BF development. DE lncRNAs show similar three expression patterns with differentially expressed protein-coding genes. Of note, cis target genes of DE lncRNAs in cluster 1 were enriched in leukocyte migration, positive regulation of locomotion, development growth, etc. These enriched GO catergories are closely linked with functional process of B lymphocytes. It is well documented that fully differentiated B lymphocytes would migrate from the BF to peripheral lymphoid organs to colonize, reproduce and perform important immune functions (5). Representative target genes of DE lncRNAs include IL16, CXCR4 and ITGA4, etc. IL-16, also known as lymphocytic chemoattractant factor (LCF), can regulate the migration and proliferation of normal leukocytes (37). The chemokine receptor CXCR4 is abundant on naive lymphocytes. It is a critical regulator of cell migration (38). ITGA4 is a member of the integrin family of transmembrane receptors, which is involved in generating signals pivotal for cell migration (39). Combining our results with the functions of these target genes, it can be speculated that lncRNAs might participate in the regulation of pigeon BF development and B lymphocyte function.

Based on our findings that immune-related mRNAs and lncRNAs were identified as differentially expressed, we speculated that, like other poultries, pigeon BF plays important roles in B lymphocyte development and differentiation. The reasons why pigeon squabs subjected to bursectomy didn’t exhibit a decline of humoral immune responsiveness were possibly related to the following facts: Firstly, unlike other poultries, parent pigeons have the capability to secrete crop milk, which contains immunoglobulin, protein and fat, etc (40, 41). By means of feeding crop milk, pigeon parents could transfer immunoglobulin to squabs (42). Secondly, emigration of B cells from bursa to periphery lymphoid tissues may start around the time of hatch (43). Thirdly, Approximately 5% of peripheral B cells in the newly hatched bird appear to be derived from extrabursal precursors (44).

Overall, in this study we presented the results of lncRNA and mRNA expression profiling in pigeon BF across three different development stages and highlighted their potential biological function. Our findings will provide insight into the transcriptional difference involving in pigeon BF development and enhance the understanding the physiological functions of squab BF.

Materials and methods

Preparation of experimental animals and tissues

White King pigeon squabs were purchased from the Feng Mao pigeon breeding plant (Mianyang, China). Thirty-nine pigeon squabs were grouped into three stages (13 replicates for each period, the age of 1 day (1d), 13 days (13d), 26 days(26d). After measuring the body weight, the pigeon squabs were anesthetized with diethyl ether and euthanatized. Then, the BF of each pigeon was isolated and weighed. BF index was calculated using the formula: BF index = BF weight/body weight ×100. And the samples for subsequent lncRNA and mRNA sequencing were frozen in liquid nitrogen and stored at -80°C refrigerator.

Histomorphological examination of BF

BFs were fixed in 4% paraformaldehyde for over 24h and dehydrated in graded ethanol. Subsequently, BFs were embedded in paraffin and subjected to a microtome (Leica, Wetzlar,Germany). Serial 5μm-thick slices were prepared and stained with H&E, and taken photos with a digital camera.

RNA extraction and high-throughput sequencing

Total RNA was respectively isolated from frozen BF samples collected from nine female pigeon squabs across the three developmental stages (three biological replicates each) by HiPure Universal RNA Kit (Magen, Guangzhou, China) according to the manufacturer’s protocol. The integrity of total RNA was evaluated using an Agilent 2100 Bioanalyzer. The RNA concentrations were assessed with a Nanodrop 2000 spectrophotometer (Thermo Scientific, Massachusetts, USA). Each total RNA sample was first treated to remove ribosomal RNA using the Ribo-Zero™ kit (Epicentre,Madison, WI, USA). The rRNA-depleted RNA was used to prepare the stranded-specific RNA-seq library following the manufacturer’s standard procedures. Subsequently, prepared cDNA libraries were sequenced by an Illumina NovaSeq 6000 platform and the paired-end reads were obtained. All the RNA-seq data have been deposited in NCBI database with accession number GSE183791.

Assembly of pigeon transcripts and identification of lncRNA

High-quality clean reads were obtained by removing reads containing adapters and low-quality reads from the raw data. These clean reads were mapped to the pigeon genome (Cliv_1.0 from NCBI) using STAR (45) (v.2.6.0c). Mapped reads were assembled to transcripts using the Cufflinks (46) (v.2.1.1) software. The transcripts (≤250 bp or FPKM < 0.1) were first excluded. The remaining transcripts were merged and compared to the reference genome by the meta-assembly tool TACO software (v0.7.3). Those transcripts unannotated by protein-coding genes (PCGs) were retained as putatively lncRNA transcripts. Subsequently, Coding Potential Calculator (CPC2) and PfamScan (47) were applied to predict the coding potential of these transcripts. After removing transcripts with coding capacities, the remaining transcripts were identified to be long non-coding RNAs.

Differential expression and cluster analysis

Quantification of mRNA and lncRNA expression in each sample was calculated in TPM (transformed transcripts per kilobase million) scale by Kallisto (48) (version 0.43.0). Next, mRNAs with TPM > 0.5 and lncRNAs with TPM > 0.1 in at least one library were considered expressed. Differential expression analysis between two stages was performed using edgeR (49)(Release 3.13) in the OmicShare tools (www.omicshare.com/tools), with |log2(fold change) |>1 and FDR <0.05. The expression patterns of differentially expressed mRNAs and lncRNAs were characterized by R package TCseq, respectively (50).

Target gene prediction of DE lncRNA

To explore the possible function of DE lncRNAs, cis target genes of lncRNAs were predicted. DE genes located 100 kb upstream and downstream of lncRNAs were extracted. Pearson’s correlation coefficients (PCC) between DE mRNAs and DE lncRNAs were calculated with the Hmisc package in R. The lncRNA-mRNA pairs with |PCC| > 0.9 and P < 0.05 were predicted as co-regulated.

Functional enrichment analysis

The DE mRNAs and cis target genes of DE lncRNAs were subjected to functional enrichment analyses by using Metascape (http://metascape.org/gp/index.html) online tool with the default parameters set (51). The GO terms and KEGG pathways with P < 0.05 were considered significantly enriched.

Protein-protein interaction network analysis of DE genes

Protein-protein interaction (PPI) networks were constructed by using the Metascape tool with default parameters. Next, MCODE algorithm was employed to identify highly connected network components from the PPI network. Functional enrichment analysis was conducted for each MCODE component, and the three best-scoring terms (by P-value) were retained as the functional description of the corresponding modules.

Quantitative real-time PCR validation

The total RNA was reverse transcribed respectively by PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Otsu, Japan) following the manufacturer’s recommendation. The qRT-PCR assays were performed on a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA) using SYBR Premix Ex Taq II (Novoprotein, Shanghai, China) by the manufacturer’s protocols. Relative mRNA and lncRNA levels were calculated by the 2^-△△Ct method. The primer sequences are in Supplementary Table 1.

Statistical analysis

BF weight and index were analyzed and compared for statistically significant difierences using one-way analysis of variance (ANOVA) in R, respectively. The significance level for the statistical analysis was P < 0.05.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI Gene Expression Omnibus (GEO) repository under GEO accession GSE183791.

Ethics statement

This study was reviewed and approved by the Institutional Animal Care and Use Committee in the College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, China.

Author contributions

Conceptualization: XW and ML. Methodology: LJ, JM, LZ, YL and QT. Validation: HG and AJ. Investigation: JW, QP and FY. Writing-original draft preparation: XW and JW. Writing-review and editing: XL and HT. Visualization: SH. Funding acquisition: ML. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the Sichuan Science and Technology Program (2021YFYZ0009), the Ya’an Science and Technology Program (21SXHZ0022).

Acknowledgments

We thank the High-Performance Computing Platform of Sichuan Agricultural University and Ya’an Big Data Industrial Park for providing computing resources and support that have contributed to these research results.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.916086/full#supplementary-material

References

1. Davison F. The importance of the avian immune system and its unique features. In: Schat KA, Kaspers B, Kaiser P, editors. Avian Immunology, 2nd ed. Boston, MA, USA: Academic Press (2014). p. 1–9.

Google Scholar

2. Liu XD, Zhang FB, Shan H, Wang SB, Chen PY. Mrna expression in different developmental stages of the chicken bursa of fabricius. Poult Sci (2016) 95(8):1787–94. doi: 10.3382/ps/pew102

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Liu XD, Zhou B, Feng XL, Cao RB, Chen PY. Bp8, a novel peptide from avian immune system, modulates b cell developments. Amino Acids (2014) 46(12):2705–13. doi: 10.1007/s00726-014-1824-x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Fellah JS, Jaffredo T, Nagy N, Dunon D. Development of the avian immune system. In: Schat KA, Kaspers B, Kaiser P, editors. Avian Immunology, 2nd ed. Boston, MA, USA: Academic Press (2014). p. 45–63.

Google Scholar

5. Zhang Y, Zhou Y, Sun G, Li K, Li Z, Su A, et al. Transcriptome profile in bursa of fabricius reveals potential mode for stress-influenced immune function in chicken stress model. BMC Genomics (2018) 19(1):918. doi: 10.1186/s12864-018-5333-2

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Liu XD, Feng XL, Zhou B, Cao RB, Li XF, Ma ZY, et al. Isolation, modulatory functions on murine b cell development and antigen-specific immune responses of Bp11, a novel peptide from the chicken bursa of fabricius. Peptides (2012) 35(1):107–13. doi: 10.1016/j.peptides.2012.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Feng XL, Liu QT, Cao RB, Zhou B, Wang FQ, Deng WL, et al. A bursal pentapeptide (Bpp-I), a novel bursal-derived peptide, exhibits antiproliferation of tumor cell and immunomodulator activity. Amino Acids (2012) 42(6):2215–22. doi: 10.1007/s00726-011-0961-8

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cooper MD, Raymond DA, Peterson RD, South MA, Good RA. The functions of the thymus system and the bursa system in the chicken. J Exp Med (1966) 123(1):75–102. doi: 10.1084/jem.123.1.75

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Nakamura K, Yuasa N, Abe H, Narita M. Effect of infectious bursal disease virus on infections produced by escherichia coli of high and low virulence in chickens. Avian Pathol (1990) 19(4):713–21. doi: 10.1080/03079459008418726

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Yamamoto Y, Glick B. A comparison of the immune response between two lines of chickens selected for differences in the weight of the bursa of fabricius. Poultry Sci (1982) 61(10):2129–32. doi: 10.3382/ps.0612129

CrossRef Full Text | Google Scholar

11. Møller AP. Successful city dwellers: A comparative study of the ecological characteristics of urban birds in the Western palearctic. Oecologia (2009) 159(4):849–58. doi: 10.1007/s00442-008-1259-8

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ciriaco E, Gagliardi ME, Cicciarello R, Germanà G, Bronzetti P. Development of the pigeon bursa of fabricius. A scanning and transmission electron microscope study. Anatomischer Anzeiger (1985) 159(1-5):55–63.

PubMed Abstract | Google Scholar

13. Coignoul F, Vindevogel H. Cellular changes in the bursa of fabricius and thymus of cyclophosphamide-treated pigeons. J Comp Pathol (1980) 90(3):395–400. doi: 10.1016/0021-9975(80)90008-0

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Thaxton JP, Grissom RE Jr. Immunological responsiveness of the pigeon (Columba livia) following neonatal bursectomy. Dev Comp Immunol (1984) 8(1):141–8. doi: 10.1016/0145-305x(84)90018-1

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Perry RB, Ulitsky I. The functions of long noncoding rnas in development and stem cells. Dev (Cambridge England) (2016) 143(21):3882–94. doi: 10.1242/dev.140962

CrossRef Full Text | Google Scholar

16. Esteller M. Non-coding rnas in human disease. Nat Rev Genet (2011) 12(12):861–74. doi: 10.1038/nrg3074

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wong NK, Huang CL, Islam R, Yip SP. Long non-coding rnas in hematological malignancies: Translating basic techniques into diagnostic and therapeutic strategies. J Hematol Oncol (2018) 11(1):131. doi: 10.1186/s13045-018-0673-6

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. String V11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47(D1):D607–D13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Stark C, Breitkreutz BJ, Reguly T, Boucher L, Breitkreutz A, Tyers M. Biogrid: A general repository for interaction datasets. Nucleic Acids Res (2006) 34(Database issue):D535–D9. doi: 10.1093/nar/gkj109

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Türei D, Korcsmáros T, Saez-Rodriguez J. Omnipath: Guidelines and gateway for literature-curated signaling pathway resources. Nat Methods (2016) 13(12):966–7. doi: 10.1038/nmeth.4077

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Li T, Wernersson R, Hansen RB, Horn H, Mercer J, Slodkowicz G, et al. A scored human protein-protein interaction network to catalyze genomic interpretation. Nat Methods (2017) 14(1):61–4. doi: 10.1038/nmeth.4083

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Glick B. The bursa of fabricius: The evolution of a discovery. Poult Sci (1994) 73(7):979–83. doi: 10.3382/ps.0730979

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wang W, Wu P, Li Y, Hou X. Genome-wide analysis and expression patterns of zf-hd transcription factors under different developmental tissues and abiotic stresses in Chinese cabbage. Mol Genet Genomics (2016) 291(3):1451–64. doi: 10.1007/s00438-015-1136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Otsubo Y, Chen N, Kajiwara E, Horiuchi H, Matsuda H, Furusawa S. Role of bursin in the development of b lymphocytes in chicken embryonic bursa of fabricius. Dev Comp Immunol (2001) 25(5-6):485–93. doi: 10.1016/s0145-305x(00)00070-7

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Nuthalapati NK, Evans JD, Taylor RL, Branton SL, Nanduri B, Pharr GT. Transcriptomic analysis of early b-cell development in the chicken embryo. Poultry Sci (2019) 98(11):5342–54. doi: 10.3382/ps/pez354

CrossRef Full Text | Google Scholar

26. Patrone L, Henson SE, Wall R, Malone CS. A conserved sequence upstream of the B29 (Ig beta, Cd79b) gene interacts with Yy1. Mol Biol Rep (2004) 31(1):1–11. doi: 10.1023/b:mole.0000013489.04734.5e

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Komatsu A, Otsuka A, Ono M. Novel regulatory regions found downstream of the rat B29/Ig-beta gene. Eur J Biochem (2002) 269(4):1227–36. doi: 10.1046/j.1432-1033.2002.02757.x

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gong S, Nussenzweig MC. Regulation of an early developmental checkpoint in the b cell pathway by ig beta. Sci (New York NY) (1996) 272(5260):411–4. doi: 10.1126/science.272.5260.411

CrossRef Full Text | Google Scholar

29. Thompson AA, Talley JA, Do HN, Kagan HL, Kunkel L, Berenson J, et al. Aberrations of the b-cell receptor B29 (Cd79b) gene in chronic lymphocytic leukemia. Blood (1997) 90(4):1387–94. doi: 10.1182/blood.V90.4.1387

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Crotty S, Johnston RJ, Schoenberger SP. Effectors and memories: Bcl-6 and blimp-1 in T and b lymphocyte differentiation. Nat Immunol (2010) 11(2):114–20. doi: 10.1038/ni.1837

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Breloer M, Fleischer B. Cd83 regulates lymphocyte maturation, activation and homeostasis. Trends Immunol (2008) 29(4):186–94. doi: 10.1016/j.it.2008.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Su H, Na N, Zhang X, Zhao Y. The biological function and significance of Cd74 in immune diseases. Inflammation Res (2017) 66(3):209–16. doi: 10.1007/s00011-016-0995-1

CrossRef Full Text | Google Scholar

33. Sadler CR, Glick B. The relationship of the size of the bursa of fabricius to antibody Production1. Poultry Sci (1962) 41(2):508–10. doi: 10.3382/ps.0410508

CrossRef Full Text | Google Scholar

34. Wu Z, Gao S, Zhao X, Chen J, Keyvanfar K, Feng X, et al. Long noncoding rnas of single hematopoietic stem and progenitor cells in healthy and dysplastic human bone marrow. Haematologica (2019) 104(5):894–906. doi: 10.3324/haematol.2018.208926

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zeng Y, Ren K, Zhu X, Zheng Z, Yi G. Long noncoding rnas: Advances in lipid metabolism. Adv Clin Chem (2018) 87:1–36. doi: 10.1016/bs.acc.2018.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Brazão TF, Johnson JS, Müller J, Heger A, Ponting CP, Tybulewicz VL. Long noncoding rnas in b-cell development and activation. Blood (2016) 128(7):e10–9. doi: 10.1182/blood-2015-11-680843

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Atanackovic D, Hildebrandt Y, Templin J, Cao Y, Keller C, Panse J, et al. Role of interleukin 16 in multiple myeloma. J Natl Cancer Institute (2012) 104(13):1005–20. doi: 10.1093/jnci/djs257

CrossRef Full Text | Google Scholar

38. Schmidt TH, Bannard O, Gray EE, Cyster JG. Cxcr4 promotes b cell egress from peyer's patches. J Exp Med (2013) 210(6):1099–107. doi: 10.1084/jem.20122574

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Hsia DA, Lim S-T, Bernard-Trifilo JA, Mitra SK, Tanaka S, den Hertog J, et al. Integrin Alpha4beta1 promotes focal adhesion kinase-independent cell motility via Alpha4 cytoplasmic domain-specific activation of c-src. Mol Cell Biol (2005) 25(21):9700–12. doi: 10.1128/MCB.25.21.9700-9712.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Engberg RM, Kaspers B, Schranner I, Kosters J, Losch U. Quantification of the immunoglobulin classes igg and iga in the young and adult pigeon (Columba livia). Avian Pathol (1992) 21(3):409–20. doi: 10.1080/03079459208418859

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Shetty S, Bharathi L, Shenoy KB, Hegde SN. Biochemical properties of pigeon milk and its effect on growth. J Comp Physiol B (1992) 162(7):632–6. doi: 10.1007/BF00296644

CrossRef Full Text | Google Scholar

42. Goudswaard J, van der Donk JA, van der Gaag I, Noordzij A. Peculiar iga transfer in the pigeon from mother to squab. Dev Comp Immunol (1979) 3(2):307–19. doi: 10.1016/s0145-305x(79)80027-0

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ratcliffe MJH, Härtle S. B cells, the bursa of fabricius and the generation of antibody repertoires. In: Schat KA, Kaspers B, Kaiser P, editors. Avian immunology, 2nd ed. Boston, MA, USA: Academic Press (2014). p. 65–89.

Google Scholar

44. Kaiser P, Balic A. The avian immune system. In: Scanes CG, editor. Sturkie's avian physiology, 6th ed. San Diego, CA, USA: Academic Press (2015). p. 403–18.

Google Scholar

45. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. Star: Ultrafast universal rna-seq aligner. Bioinf (Oxford England) (2013) 29(1):15–21. doi: 10.1093/bioinformatics/bts635

CrossRef Full Text | Google Scholar

46. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of rna-seq experiments with tophat and cufflinks. Nat Protoc (2012) 7(3):562–78. doi: 10.1038/nprot.2012.016

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The pfam protein families database. Nucleic Acids Res (2012) 40(Database issue):D290–301. doi: 10.1093/nar/gkr1065

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic rna-seq quantification. Nat Biotechnol (2016) 34(5):525–7. doi: 10.1038/nbt.3519

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Robinson MD, McCarthy DJ, Smyth GK. Edger: A bioconductor package for differential expression analysis of digital gene expression data. Bioinf (Oxford England) (2010) 26(1):139–40. doi: 10.1093/bioinformatics/btp616

CrossRef Full Text | Google Scholar

50. Jun M, Gu L. Tcseq: Time course sequencing data analysis (2021). Available at: https://bioconductor.org/packages/devel/bioc/vignettes/TCseq/inst/doc/TCseq.pdf.

Google Scholar

51. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun (2019) 10(1):1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar