Research Note: Circular RNA expressing in different developmental stages of the chicken bursa of Fabricius.

The bursa of Fabricius plays an essential role in B lymphocyte development, which is controlled not only by proteins but by noncoding RNA. Circular RNA (circRNA) are expressed in diverse tissues in eukaryotes. To acquire a deeper perception of the molecular mechanism of bursal development, RNA sequencing was used to identify the circRNA during varied evolving stages of the chicken bursa of Fabricius. We identified 13,689 circRNA. All these circRNA were originated from 4565 chicken genes. Among them, only 1 circRNA was yielded from those 4131 parental genes, and 2 or more circular isoforms were generated from the remaining genes. There were 27 circRNA found to be differentially expressed between the embryonic day 20 and day 2 developmental stages. The 5 isoforms of immunoglobulin lambda-like polypeptide 1 circRNA were tested to validate the RNA sequencing data, and their targeted genes were also analyzed with quantitative reverse transcription PCR. These data indicate that cirRNA are abundant and essential during bursal development and may play essential roles in the development of bursa of Fabricius.

Introduction

B cells of the chicken develop in the bursa of Fabricius (BF), a special organ for birds (Davison et al., 2008). As the BF is essential for early B lymphocyte development (Glick et al., 1956, Cooper et al., 1966, Lydyard et al., 1976), B cell is named after BF. Hence, it is the BF that provides a crucial model for studies on basic immunology.

Three distinct temporal stages distinguished during chicken B cell development (Davison et al., 2008). First, the prebursal stage in which hematopoietic precursor cells commit to the B cell lineage and spread to other organs (Ratcliffe et al., 2005). Second, it is the bursal stage. The precursor cells migrate between embryo day (ED) 8 and ED 12 continuously (Douarin et al., 1975). Third, the B cell antigen receptor (BCR) is diversified by gene conversion from ED 15. Immature B cells, combined with the matured BCR, begin to move to the periphery lymphoid organs around hatch (Paramithiotis and Ratcliffe, 2010).

Noncoding RNA exist in varied cells and predominate multifarious biological functions (Sarah et al., 2012). Circular RNA (circRNA) were defined as a new class of RNA till recent progress in RNA-sequencing (RNA-Seq) technology. So far, they are widely recognized in eukaryotic tissues and cells in human beings, plants, and mice. (Julia et al., 2013, Sebastian et al., 2013). It is still largely unknown concerning the biological functions of circRNA, but recent studies indicate that some circRNA are essential in many physiological and pathologic conditions through regulation of gene expression at multiple levels (Li et al., 2018).

Here, RNA-Seq is used to identify circRNA during chicken BF development to investigate the functions of circRNA in the chicken BF.

Material and methods

Ethics Committee

As for the chicken BF in this work, the experimental protocols were approved by the Animal Care and Use Committee of Nanjing Agricultural University, China (Approval No. 200709005).

Sample Preparation

Fertilized specific pathogen–free eggs of White Leghorn chickens were supplied by QYH Biotech Company (Nanjing, China). It was at QYH Biotech Company (Nanjing, China) that specific pathogen–free chickens were hatched and reared. Organ samples were picked up from ED 20 (ED20) (n = 3) and from birds at day 2 (D2) (n = 3) after being hatched. Each treatment was run in 3 replicates.

Whole-Transcriptome Sequencing and Data Analysis

The RNA which was extracted using TRIzol reagent (Invitrogen) with RIN >6.0 was applied to compose rRNA depletion library (NEBNext Ultra Directional RNA Library Prep Kit) in accordance with the manufacturer's directions. The whole-transcriptome sequencing data were obtained with a HiSeq Sequencer. The data were filtered (removing the adapter sequences, reads with >5% ambiguous bases (noted as N) and low-quality reads containing more than 20% of bases with a quality score llower than 20. The reads were then mapped to the chicken genome (GRCg6a NCBI) using the HISAT2 alignment program (Daehwan et al., 2015). The HTSeq program was used to determine the gene count of mRNA and long non-coding RNA (Simon et al., 2015).

Circular RNA Prediction and Data Analysis

Circular RNA was previewed on the basis of the sequencing data with ACFS circRNA prediction pipeline (Xintian et al., 2015). Unmapped reads were collected by using BWA-MEM (bwa mem -t 1 -k 16 -T 20) for circRNA identification: the junction of head-to-tail was recognized. The maximum score concerning the RNA splice site strength was counted using MaxEntScan33, the filtering standard of which is larger than or equal to 10. Reads that mapped to the circRNA back-splicing junction (with an overhang of at least 6 nt) were counted.

Differential Expression Analysis and Competing Endogenous RNA Relation Prediction

The differential expression analysis of the RNA and circRNA was conducted with the DESeq software (Anders and Huber, 2010) with the criteria: Fold Change >2; FDR <0.05. Miranda package (Enright et al., 2003) was introduced to make the prediction of microRNA (miRNA) target on 3′ untranslated region of disparately expressed mRNA and the full length sequence of circRNA, which was disparately expressed. Competing endogenous RNA relations was shaped by a pair of circRNA and mRNA influenced by the same miRNA members.

Co-expression Analysis

We present gene co-expression networks to find the relations among genes and circRNA. Gene co-expression networks were constructed in accordance with the normalized expression values of genes which were chosen from those in consequential pathway terms and GO terms and differentially expressed circRNA (Miguel Angel et al., 2007). For each pair of genes, the Pearson correlation was calculated, and the significant correlation pairs were selected under a false discovery rate (FDR) < 0.05 to shape the co-expression network. Moreover, to explore the network property, k-cores in graph theory were imported in this work as a method to simplify graph topology analysis, indicating the central status of a circRNA in genes and circRNA (Ravasz et al., 2002).

Functional Analysis

We used gene ontogeny analysis and pathway analysis with the GO database (download from NCBI, UniProt, and AmiGO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) database, respectively, for functional annotation. The Fisher’s exact test was used to define the significance P-value (Draghici et al., 2007) and BH test was used to calculate the FDR (Benjamini and Hochberg, 1995).

Real Time PCR

For validation of gene expression, total RNA was prepared from bursas collected at ED18, ED20, and D2 with the TRIzol Reagent (Invitrogen). Real-time PCR was displayed in a LightCycler instrument (Eppendorf, German) in a total volume of 20 μL with the application of the SYBR Premix Ex Taq (Perfect Real Time) kit (TaKaRa Bio Inc, Shiga, Japan) in accordance with the manufacturer's directions. The Supplementary Table 1 shows the primer sets from selected genes.

Results and discussion

Overall Quality Parameters of the Circular RNA-Seq Data

Various researches have indicated that circRNA are widely distributed and are found in many eukaryotes, with mammals in particular. Nevertheless, little is known concerning circRNA in the BF of chickens. In this work, the genome-wide recognition was performed, and the function analysis of circRNA in the BF was conducted accordingly.

To detect quantitative differences of circRNA during chicken bursa development, the BF at ED20 and D2 representing the significant developmental phases of bursal B cell was used.

Illumina high-throughput second-generation sequencing produced 71746510, 79503310, and 81081288 and 65604606, 80904166, and 69846244 clean reads for the ED20 and D2 stages, respectively (Supplementary Table 2). Approximately 92–93.3% of those reads (Supplementary Table 3) were mapped equivalently for all samples to the chicken genome. From the six bursa tissues, a total of 13,689 circRNA were detected in 2 stages (Supplementary Table 4). The detected circRNA were distributed in each chromosome, except chromosome MT (Figure 1A). All these circRNA were derived from 4565 chicken genes. Among them, only 1 circRNA was produced by 4131 parental genes, and 2 or more circular isoforms were yielded by the remaining genes (Figure 1B).

Figure 1

Annotation of chicken bursa circRNA. (A) Distribution of sequencing reads of circRNA on chicken chromosome. (B) Distribution of circRNA among genes. Abbreviations: circRNA, circular RNA; ED, embryonic day.

Identification of Differentially Expressed circRNA Among the Different Development Stages

To acquire a deeper perception of the biological mechanism of the chicken bursa development, it is critical to identify the disparate expression circRNA (DEcircRNA) during different phases. There were 27 DEcircRNA that were detected between the ED20 and D2 groups when FDR 0.05 was used as cutoff values (Supplementary Table 5). Of these DEcircRNA, 10 were upregulated, whereas 17 were downregulated in ED20 compared with those in D2. Figure 2A shows the expression heat map on the matrix of variance-stabilized data.

Figure 2

Circular RNA expression. (A) Expression heat map on the matrix of variance-stabilized data for overall circRNA expression (green, no or very low expression; red, high expression). (B) KEGG pathway enrichments of ED20 vs. D2. (C) GO biological process terms of ED20 vs D2. (D) GO molecular functional terms of ED20 vs. D2. Abbreviations: circRNA, circular RNA; GO, gene ontology; ED, embryonic d; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Gene ontology analysis is a functional analysis which associates differentially expressed genes with gene ontology categories (Ashburner et al., 2000). The biological process categories analysis showed that parental gene functions were involved in regulation of membrane depolarization during action potential, regulation of DNA-templated transcription, initiation, “de novo” GDP-L-fucose biosynthetic process, cytoplasm organization, negative regulation of intrinsic apoptotic signaling pathway by p53 class mediator, and nucleolus organization (Figure 2C). The molecular functional categories analysis showed that parental gene functions were involved in nicotinamide adenine dinucleotide phosphate + binding, microfibril binding, GDP-mannose 4,6-dehydratase activity, hedgehog family protein binding, oxidoreductase activity, protein phosphatase type 4 regulator activity, and phosphoprotein phosphatase activity (Figure 2D).

Then, we mapped the differentially expressed genes to Kyoto Encyclopedia of Genes and Genomes terms to discover those genes involved in metabolic or signal transduction pathways that were significantly enriched (Figure 2B). The pathway analysis revealed that parental gene functions were involved in fructose and mannose metabolism, amino sugar and nucleotide sugar metabolism, Hedgehog signaling pathway, cardiac muscle contraction, gap junction, gonadotropin-releasing hormone signaling pathway, transforming growth factor-beta signaling pathway, vascular smooth muscle contraction, Rap1 signaling pathway, Ras signaling pathway, calcium signaling pathway, cytokine–cytokine receptor interaction, Regulation of actin cytoskeleton, focal adhesion, MAPK signaling pathway.

Putative Functions of Chicken circRNA as miRNA Sponges

Many reports showed that circRNA can function as miRNA sponges or competitive (Westholm et al., 2014, Fang et al., 2015). We identified the miRNA-binding capabilities of our differentially expressed circRNA and matched them to their potential target miRNA (top 5 showed in Supplementary Table 6). These data were then applied to form an interact network map (Figure 3). It was observed that the 5 isoforms of immunoglobulin lambda-like polypeptide 1 (IGLL1) circRNA showed the most interactions with target miRNA. These DEcircRNA could be regarded as candidate circRNA for further research in bursal development.

Figure 3

Interaction network assigned with the different expression circRNA and their potential target miRNA. The top 10 differentially expressed circRNA and their potential target miRNA were used for their interaction network analysis. Abbreviations: circRNA, circular RNA; miRNA, microRNA.

The miR-34a-5p was found to be essential for megakaryopoiesis and the production of monocytes (Bianchi et al., 2017). Interestingly, it was found that circRNA was derived from the IGLL1 gene–harbored binding sites for miR-34a-5p. A total of 5 isoforms of IGLL1 circRNA were identified and were differentially expressed in the various stages of bursal development (Figure 4A and Supplementary Table 7). The interaction of circLGLL1 with miR-34a-5p will be the subject of future studies.

Figure 4

The biological validation of differential expression. Heatmap of IGLL1 circRNA and B cell differentiation–related genes were shown. The data were the fold induction compared with the lower expressed group. (B) The biological validation of differential circRNA expression. (C) Expression of PAX5 and IGLL1genes. qRT-PCR reactions were run in triplicate and presented as means ± SD, ∗∗P < 0.01. Abbreviations: circRNA, circular RNA; IGLL1, immunoglobulin lambda-like polypeptide 1; PAX5, paired box 5; qRT-PCR, quantitative reverse transcription PCR.

Validation of Chicken circRNA

Biological validation was conducted with the use of quantitative reverse transcription PCR. All tested cirRNA demonstrated the same tendency with the original RNA-Seq data set (Figure 4B).

For IGLL1 and paired box 5 (PAX5), three developmental stages (ED18, ED20, and D2) were selected from an independent set of animals to validate the RNA-Seq data. The expression patterns of these genes were consistent with the original RNA-Seq data set (Figure 4C).

During the early development of B cell, cells that are at the pre-BI stage express IGLL1, the gene encoding λ5, serving as part of the pre-B cell receptor complex (pre-BCR). In pre-BII stage, cell downregulated pre-BCR expression (Thompson et al., 2007). The IGLL1 gene expression is downregulated at the pre-Bcell stage in mice (Alexander et al., 2008). In our results, we found that IG11 was also downregulated (Figures 4A, 4C). The amounts of Ikaros and EBF expression remained the same, consistent with the data previously reported (Sabbattini et al., 2014). For cirRNA, we found that expression of all 5 isoforms of IGLL1 circRNA were downregulated from the ED20 stage to D2 stage of bursal development. The decrease of these circRNA may decrease the competition with mRNA for miRNA binding and thus downregulate the expression of IGLL1.

In conclusion, we conducted genome-wide identification of chicken circRNA from the BF with RNA sequencing method and found they are substantial and disparately expressed in various chicken bursa development. Most of the differentially expressed circRNA harbored miRNA-binding sites and could play essential roles in bursal development.