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

The temporal expression patterns of brain transcriptome during chicken development and ageing

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

Abstract

Background

The transcriptional profiles of mammals during brain development and ageing have been characterized. However the global expression patterns of transcriptome in the chicken brain have not been explored. Here, we systematically investigated the temporal expression profiles of lncRNAs and mRNAs across 8 stages (including 3 embryonic stages, 2 growth stages and 3 adult stages) in the female chicken cerebrum.

Results

We identified 39,907 putative lncRNAs and 14,558 mRNAs, investigated the temporal expression patterns by tracking a set of age-dependent genes and predicted potential biological functions of lncRNAs based on co-expression network. The results showed that genes with functions in development, synapses and axons exhibited a progressive decay; genes related to immune response were up-regulated with age.

Conclusions

These results may reflect changes in the regulation of transcriptional networks and provide non-coding RNA gene candidates for further studies and would contribute to a comprehensive understanding of the molecular mechanisms of chicken development and may provide insights or deeper understanding regarding the regulatory mechanisms of age-dependent protein coding and non-protein coding genes in chicken. In addition, as the chicken is an important model organism bridging the evolutionary gap between mammals and other vertebrates, these high resolution data may provide a novel evidence to improve our comprehensive understanding of the brain transcriptome during vertebrate evolution.

Background

Large-scale genomic studies have revealed that the majority of cellular RNAs are transcribed as non-protein coding RNAs (ncRNAs) [1,2,3,4]. Long non-coding RNAs (lncRNAs) are a subgroup of RNAs longer than 200 nucleotides (nt); many of them are 5′ capped, alternatively spliced and polyadenylated as mRNAs are [5]. LncRNAs are newly recognized regulating molecules and can be classified as intergenic lncRNAs (lincRNAs), intronic lncRNAs, antisense lncRNAs or sense-overlapping lncRNAs based on their genomic location [6]. Albeit mRNA-like lncRNAs have limited protein-coding potential, they function in many biological processes and play roles as signaling molecules, scaffolds, guides and decoys [7]. Transcriptional dynamics has been suggested to be a major contributor for brain structure development, as well as to the ageing [8]. Most lncRNAs are less conserved and expressed at significantly lower levels, as well as in more restricted temporal expression patterns, than mRNAs are, indicating that lncRNAs are of considerable importance in modulating various biological functions, especially development and cellular differentiation [9,10,11].

A broad range of transcriptional profiles in fish [9, 12], mouse [10, 13,14,15,16], monkeys [17,18,19] and humans [1, 20,21,22,23,24,25,26] during brain development and ageing have been reported in recent years. It has been demonstrated that the encoded information underlying development lies in regulatory elements that define different gene expression programs, which means the regulatory factors scale quadratically with genome size, the untranslated regions (UTRs) in mRNAs increase with the developmental complexity of organisms, and the non-protein coding intronic and intergenic transcripts also contain an expanded regulatory framework to control gene expression during differentiation and development [20]. Although remarkable achievements have been made in characterizing the transcriptional underpinnings of brain development and ageing in vertebrates, the precise mechanisms of lncRNAs in defining the complexity of brain functions remain unclear, particularly in birds. Although the chicken is now one of the most versatile experimental systems available, as an established model organism for the study of vertebrate development, the chicken is under evaluated in the study of developmental transcriptomics [27].

To systematically investigate the expression profiles of lncRNA and mRNA transcriptome in different stages, we performed high-throughput RNA sequencing to interrogate their temporal expression patterns in the chicken cerebrum. The present study provided an initial lncRNA catalogue of the chicken cerebrum, which improves our understanding of the temporal expression profile during development and ageing, and may provide insights or deeper understanding regarding the regulatory mechanisms of age-dependent genes in chicken. In addition, this lncRNA catalogue helps bridge the evolutionary gap between mammals and other vertebrates.

Results

Global identification of lncRNAs in the chicken cerebrum

To systematically investigate lncRNAs and their temporal expression profiles during development and ageing in the chicken cerebrum, a total of 21 cDNA libraries from 21 female Tibetan chickens at 8 life stages were constructed. The sampling times broadly spanned the life of the farm-raised chicken, representing 3 embryonic periods (embryonic days 12, 16, and 20: E12, E16 and E20), 2 growth periods (100 and 300 days of age: D100 and D300), and 3 adult periods (from young adults to ageing adults--1, 3, and 5 years of age: Y1, Y3 and Y5). Illumina sequencing yielded a total of 2,675,315,228 raw reads and 406.56 Gb of raw data with 150-bp paired-end RNA-Seq, 2,615,649,932 clean reads were obtained after quality control, corresponding to an average of 18.91 Gb of high-quality data per sample. Among them, 91.87–94.61% were mapped to the chicken genome (Gallus_gallus-5.0) through HISAT2 2.1.0 [28] (Additional file 1 Table S1).

Genomic characterization and classification of lncRNAs

To compare the genomic characteristics of lncRNAs and mRNAs, we identified a total of 39,907 putative lncRNAs and 14,558 mRNAs that were expressed in more than one biological replicate (FPKM > 0). Subsequently, the genome distribution of lncRNAs and mRNAs was investigated, of which there were 36,427 (86.04%) transcripts and 12,254 (93.02%) mRNAs were assembled on chromosomes (Fig. 1a and Additional file 1 Table S2). Most of the lncRNAs and mRNAs of chicken tended to distribute in large chromosomes (No. 1–10 and Z chromosome), while a small proportion of transcripts distributed in micro chromosomes. In addition, the number of lncRNAs and mRNAs were highly correlated with chromosome size (r = 0.9850, 0.9677, respectively). Illustrating that the number of transcripts distributed on chromosomes was proportional to the size of chromosome. Based on further comparisons, we found that chicken lncRNAs are approximately 1/2 the length of mRNAs (mean length of 1122 nt for lncRNAs and 2029 nt for mRNAs); moreover, lncRNAs have fewer exons and are expressed at lower levels than mRNAs (Fig. 1b-d).

Fig. 1
figure 1

Genomic characterization of lncRNAs. a Chromosome distribution of mRNA and lncRNA transcripts identified in chicken cerebrum. The blue line represents the size of chromosome (the ratio of each chromosomal size to the total genomic size). b The distribution of transcript length of lncRNAs and mRNAs. c Exon number distribution of lncRNAs and mRNAs. d Comparison of the expression levels of lncRNAs and mRNAs. The lines of the whiskers in the box represent median lines. e The classification of lncRNAs

Full size image

We identified 36,663 lncRNAs and classified them by genomic location; 54.24% of them were classified as intergenic lncRNAs not overlapping with any genes, 16.27% as sense exon overlapping lncRNAs, 8.62% as antisense exon overlapping lncRNAs, 6.80% as sense intron overlapping lncRNAs and 14.08% as antisense intron overlapping lncRNAs (Fig. 1e and Additional file 1 Figure S1). Similarly, the subgroups of lncRNAs (intergenic, exonic and intronic lncRNA) distributed on chromosomes in proportion to the chromosome size (Additional file 1 Table S2 and Figure S2). Furthermore, we performed functional enrichment analysis on both exon overlapping lncRNAs and intron overlapping lncRNAs, terms related to neuron, axon and synapse were found to be enriched in these genes (Additional file 2). For instance, GO: 0043524 (negative regulation of neuron apoptotic process), GO: 0048666 (neuron development) and GO: 0007269 (neurotransmitter secretion). Axon activities including GO: 0007411 (axon guidance), GO: 0008089 (anterograde axonal transport), GO: 0031290 (retinal ganglion cell axon guidance), GO: 0019896 (axonal transport of mitochondrion), GO: 1904115 (axon cytoplasm) and GO: 0071679 (commissural neuron axon guidance). Terms involved in synaptic regulation contained GO: 0045211 (postsynaptic membrane), GO: 0051965 (positive regulation of synapse assembly), GO: 0007416 (synapse assembly), GO: 0048813 (dendrite morphogenesis). KEGG pathways were primarily involved in cell proliferation and differentiation, development such as “Insulin signaling pathway”, “MAPK signaling pathway”, “ErbB signaling pathway”, “Hedgehog signaling pathway”; “Neuroactive ligand-receptor interaction” and “Calcium signaling pathway” were related to neuron. These results indicated that the intron and exon overlapping lncRNA functions were consistent with the biological functional activities of the genes enriched in the cerebrum.

Expression profiles of mRNAs and lncRNAs

Our RNA-Seq experiments included a broad time series of sampling which allowed us to follow the expression dynamics of mRNAs and lncRNAs as development proceeded. Based on the expression dynamics of mRNAs and lncRNAs, the Pearson correlation between each pair of samples was calculated. The results showed that 21 samples were grouped into three broad clusters of mRNA expression profiles (Fig. 2a): the first cluster (green) represented the embryonic stages (E12, E16 and E20), the second cluster (blue) meant the growth stages (D100 and D300) and the last cluster (red) was the adult stages (Y1, Y3 and Y5). Chickens included in the red cluster had reached a growth plateau and even begun the ageing process. While the expression profile of lncRNAs was not as markedly varied as that of mRNAs, the embryonic and growth stages were clustered, and the adult stages were congregated (Fig. 2c). In addition, the correlation power of mRNAs at consecutive time points exhibited strong time-dependent clustering. Meanwhile, the correlation power of lncRNAs between stages of lncRNAs was much weaker (Fig. 2b, d). Subsequently, we calculated the Shannon entropy (H) value [29] to further test the hypothesis; this parameter can be used to measure the specificity of gene expression during developmental stages. Compared with mRNAs, namely, lincRNAs, intron overlapping and exon overlapping lncRNAs, exhibited increased temporal specificity (Fig. 2e).

Fig. 2
figure 2

The expression profiles of mRNAs and lncRNAs. a, c The expression profiles of mRNAs (upper panel) and lncRNAs (lower panel). The heat map shows the expression profile of mRNAs and lncRNAs. The top and left panel is the sample and gene tree constructed by Pearson Correlation; the value represents the log2 transformed values of (FPKM+ 1). b, d Dynamic changes in expression profile of mRNAs and lncRNAs. The value represents the pairwise Pearson Correlation. The correlation between every two samples was calculated by log2-transformed (FPKM+ 1) gene expression values. e The distributions of Shannon entropy-based temporal specificity scores that were calculated for distinct classes of mRNAs and lncRNAs. f Age-dependent clusters with significantly enriched GO terms and pathways map. Hierarchical clustering analysis grouped the age-dependent genes into 7 clusters. g Temporal profiles of age-dependent DEGs (upper panel) and DE lncRNAs (lower panel). Modules in different colors represent different temporal expression patterns which were significantly enriched by STEM analysis (non-significant modules were not shown). Numbers in the top indicates the module number and the lower numbers indicate numbers of DEGs and DE lncRNAs in each module

Full size image

Age-dependent transcriptome and gene expression analysis

To evaluate differentially expressed lncRNAs and mRNAs (DE lncRNAs and DEGs) at irregularly spaced sampling times, we compared the expression levels at adjacent time points with the threshold of FC(Fold Change) ≥ 2 (or ≤ 0.5) plus a Bonferroni-adjusted P-value ≤ 0.05. Overall, we compared 7 pairs and identified 1892 unique lncRNAs and 4499 unique mRNAs as significantly differentially expressed during development and ageing. In the lncRNAs comparison libraries, 1122 and 1059 lncRNAs were found to be up- and down-regulated, respectively; a total of 2985 up-regulated and 2993 down-regulated DEGs were identified in mRNAs comparison groups during chicken cerebrum development and ageing. The greatest number of DE lncRNAs between adjacent time points was observed at Y5 vs Y3, followed by Y1 vs D300 and E16 vs E12, while the most DEGs were detected at Y1 vs D300, followed by E20 vs E16 and E16 vs E12 (Additional file 1 Figure S3).

To further investigate the biological functional variations of these DEGs in development and ageing, we used 4499 DEGs to perform hierarchical clustering analysis. The sample tree was clustered into 3 subgroups by age, and the gene tree was partitioned into 7 clusters, which were labeled with different colors (Fig. 2f). Subsequently, we performed functional enrichment analysis based on the DEGs of the 7 clusters. Generally, DEGs expressed highly in embryonic periods were enriched in the terms “cell division”, “insulin secretion” and “axon guidance”; the most relevant KEGG pathways were “glutathione metabolism”, “Notch signaling pathway” and “cell cycle” (the first three clusters). DEGs in growth periods were enriched in the terms “MAPK signaling pathway”, “GnRH signaling pathway” and “Calcium signaling pathway” (the two middle clusters). Adult periods were mainly enriched in immune response, phagosome and defense response to virus (the last 2 clusters) (Additional file 3).

Temporal profiles of gene expression

We used STEM to group DEGs with similar temporal profiles and found that 4340 mRNAs and 1218 lncRNAs were classified into time-series expressed patterns. For mRNAs, a total of 1922 (44.29%) genes were partitioned into 9 significantly enriched temporal profiles (P ≤ 0.05). Two models exhibited monotonic trajectories: the pink and blue-violet models showed a gradual increase and decrease with age, respectively. Meanwhile, a portion of genes were assigned to models that reversed their trajectory over time: the yellow model presented a U-shaped trajectory with a minimum at D100 to D300 and the orange model showed a Bell-shaped trajectory with a peak at D300. The remaining models had more complex temporal profiles, some with an initial increase (green, light purple and jade models) and some with an initial decrease in their expression (mauve and light brown models) (Fig. 2g). As for lncRNAs, 206 (16.91%) transcripts were assigned to significantly enriched model profiles (P ≤ 0.05) (Fig. 2g). To investigate the biological functions of these significantly enriched mRNAs, we performed functional enrichment analysis on aforementioned genes in 9 models. Interestingly, the terms “cytokine receptor activity” (GO: 0004896), “phagosome” and “ECM-receptor interaction” were enriched in increasing genes; “cell cycle”, “MAPK signaling pathway”, “p53 signaling pathway”, “mitotic nuclear division” (GO: 0007067), “axon guidance” (GO: 0007411), “axon hillock” (GO: 0043203) and “cell proliferation” (GO: 0008283) were enriched in decreasing genes. In addition, “activation of GTPase activity” (GO: 0090630) and “regulation of axonogenesis” (GO: 0060070) were enriched in DEGs from Bell-shaped temporal profiles; “canonical Wnt signaling pathway” (GO: 0060070) and “growth” (GO: 0040007) were enriched in DEGs from U-shaped temporal profiles (Additional file 4).

Construction of a co-expression network

Predicting the putative function of lncRNAs through sequence characteristics is still a challenge because of the lack of annotated features. WGCNA can be applied to associate lncRNAs with functionally annotated mRNAs because genes and transcripts with similar expression patterns are often related in biological functions and can be attributed to the same module. We first used DESeq2 to screen out mRNAs and lncRNAs; secondly, we conducted WGCNA to explore the function of lncRNAs, a total of 11,200 mRNAs and 955 lncRNAs were selected by screening to construct the co-expression network. Ultimately, 12 modules were identified; the top 2 modules accounted for 74.02 and 70.26% of the total mRNAs and lncRNAs, respectively. We therefore considered a total of 8323 mRNAs were targeting genes to 671 lncRNAs. Furthermore, the co-expressed genes of the top 2 modules (8323 mRNAs) were submitted to DAVID; the functional categories showed that these co-expressed genes were significantly enriched in the terms “myelin sheath” (GO: 0043209), “dendrite morphogenesis” (GO: 0048813) and “axonogenesis” (GO: 0007409); as well as nuclear components such as “nucleolus” (GO: 0005730), “nucleoplasm” (GO: 0005654), “chromatin organization” (GO: 0006325). Consistently, KEGG analysis showed an extensive enrichment of co-expressed genes involved in pathways such as “spliceosome”, “ribosome”, “cell cycle”, “insulin signaling pathway” and “mTOR signaling pathway”. It is worth noting that, pathways related to defense, detoxication and antiageing were significantly enriched in adult stages: “positive regulation of telomere maintenance via telomerase” (GO: 0032212). “Toxin transport” (GO: 1901998), “damaged DNA binding” (GO: 0003684), “wound healing” (GO: 0042060), “glutathione metabolic process” (GO: 0006749) and “response to UV” (GO: 0009411) (Additional file 5).

Validation of the gene expression by qPCR

To validate the sequencing data, we selected 10 mRNAs and lncRNAs representing different temporal profiles to examine the expression patterns at each developmental stage. Three DEGs (TMSB15B, BASP1, GAP43) were decreasing pattern, 2 DEGs were increasing (FGF13, CACNA2D1), Bell-shaped (CAMK2A, GRIN1) and U-shaped (TTR, TMSB15B) patterns, respectively, and 1 DEG (PLP1) came from the initial-up expressing pattern. Correspondingly, 3 DE lncRNAs in pink module (TU56811, TU178751, TU167089), 3 in green module (TU99158, TU26292, TU50759), and 4 in blue module (TU38242, TU45029, TU166039, TU6367) were selected from the temporal expression profiles of lncRNAs in Fig. 2g. The relative expression of 10 mRNAs and lncRNAs detected by qPCR were compared with the transformed log2 (FPKM+ 1) values of RNA-Seq. The results showed that the altered expression patterns of 10 mRNAs and lncRNAs at each stage were consistent with the RNA-Seq data (Fig. 3), illustrating the reliability of our RNA-Seq data and guaranteed the accuracy of stringent pipeline was used to identify the transcripts.

Fig. 3
figure 3

Validation of 10 expressed DEGs and DE lncRNAs by qPCR. The X-axis represented 8 developmental stages. The Y-axis indicated the relative expression of each gene; blue lines were log2 (FPKM+ 1) values of RNA-Seq and the red lines were relative expression of qPCR

Full size image

Discussion

Genomic characterization of lncRNAs in chicken cerebrum

Although the brain lncRNA profile of human, macaque, mice and fish have been well characterized in recent years, the function of lncRNAs in chicken brain development and ageing remained unclear. Here, we systematically investigated the lncRNA profile of the cerebrum in female Tibetan chicken by RNA-Seq spanning 8 stages from embryonic (E12) to ageing (Y5) individuals. We obtained 406.56 Gb of raw data from 21 libraries in total; the chromosomal distribution of transcripts illustrated that the number of transcripts distributed on chromosomes was proportional to the size of chromosome; not surprisingly, several key genomic characteristics of the chicken cerebrum that lncRNAs were shorter, less conserved and expressed at lower levels than mRNAs were consistent with results in zebrafish and mammals [2, 8, 9, 17, 18]. Furthermore, we identified approximate 2000 (54.24%) lincRNAs (long intergenic lncRNAs) in our study, accounting for more than half of the total chicken lncRNAs, which were consistent with previous studies [30, 31]. Functional enrichment analysis indicated that the functions of intron and exon overlapping lncRNAs were consistent with the biological functional activity of the genes that were enriched in the cerebrum, and that brain transcriptomes are conserved across vertebrates.

Atlas of gene expression in chicken cerebrum

The expression profile of mRNAs of the chicken cerebrum was age-dependent and lncRNA has a more restricted temporal expression than mRNA, suggesting that there might be distinctive patterns of age-related functional changes across cerebrum developmental stages. Further findings on age-dependent transcripts showed that DEGs were grouped into 7 sub clusters. The first 3 clusters mainly highly expressed in embryonic stages, involving in the formation of neuronal network and development. For example, HDAC2 encodes histone deacetylase which can form transcriptional repressor complexes by associating with many different proteins, and play an important role in transcriptional regulation, cell cycle progression and developmental events. Previous reports have indicated that HDAC2 is critical for oocyte development and reciprocally, it also can maintain normal chromatin structures and has a unique role by controlling the fate of neural progenitors during normal brain development [32, 33]. Activity of protein product encoded by PAX6 is essential in the development of neural tissues, particularly the eye. Pax6 can establish distinct ventral progenitor cell populations and control the identity of motor neurons and ventral interneurons [34, 35]. The 2 middle clusters were up-regulated at growth stages. It was noted that, GnRH signaling pathway were enriched in this stage alone, the GnRH secreted from the hypothalamus acts upon its receptor in the anterior pituitary to regulate the production and release of the gonadotropins, LH and FSH [36, 37], which stimulated the hormones secretion and egg formation in ovary. This result was in accordance with the egg laying traits in Tibetan chicken that they began to lay eggs at the age of 240 to 270 days and reached egg laying peak at the age of 300 days [38]. Besides, GnRH could activate the signaling of PCK (protein kinase C) and lead to transactivation of the epidermal growth factor (EGF) receptor and activation of mitogen-activated protein kinases (MAPKs), resulting in activation of transcription factors and rapid induction of early genes [37]. This coincided with the fact that MAPK signaling pathway was significantly enriched in growth stage. Mainly up-regulating genes in the last 2 clusters of adult stages exerted functions in defense mechanisms and immune responses. CX3CL1 is broadly expressed in brain, it plays a role in a wide range of diseases, including cancer, vasculitis, neuropathies, atherosclerosis, inflammatory diseases, and in human immunodeficiency virus infections [39, 40]. The protein encoded by MARCO is part of the innate antimicrobial immune system [41], and modulates inflammatory responses against virus infection [42]. All these findings elucidated that biological functions of highly expressed genes in different clusters accorded with the development and growth rules.

Monotonic temporal expression patterns

The progressively increasing gene expression profile involved in phagosome and binding, which meant that organisms were undergoing neurodegeneration, anti-inflammation and innate immune mechanisms were activated as the neurons aged [12, 43, 44]. Meanwhile, the temporal profile decreased with time, of particular interest, were the down-regulation of cell proliferation and mitotic nuclear division and the degeneration of neuronal projections and axonal activity. In addition, chromatin and extracellular matrix remodeling occurred in this gene cluster. Similarly, neuronal functional deterioration and synaptic down-regulation were differentially represented in the ageing primate brain [45]. All the information indicated that a progressive accumulation of immune function and neuronal deterioration were conserved hallmarks of brain ageing.

Temporal inversions

Interestingly, a substantial proportion of DEGs reversed their expression trends at D300 and exhibited U-shaped or Bell-shaped temporal profiles. The U-shaped model applied to the “canonical Wnt signaling pathway” and “glutathione metabolism”; Wnt signaling is involved in organogenesis, and the inhibition of Wnt signaling is an essential factor in the late stage of body profile formation in vertebrates [46]. Glutathione is helpful in maintaining a normal immune system and has an antioxidant effect [2], hinting that the Wnt signaling pathway might be mapped to genes down-regulated from E12 to D100, while glutathione metabolism probably corresponds to genes up-regulated from D300 to Y5. The Bell-shaped expression model exhibited a peak at D300 and a sub-peak at E16; this gene cluster was enriched for “GTPase activity”, “regulation of axonogenesis” and “neuronal synaptic plasticity”. GTPase participated in protein biosynthesis and translocation, cell cycle control and cell differentiation [47]. Axonogenesis is a subfield of neural development concerning the process by which neurons send out axons to reach the correct targets [48]. In neuroscience, synaptic plasticity is the ability of synapses to strengthen or weaken over time, in response to increases or decreases in their activity; this phenomenon is one of the important neurochemical foundations of learning and memory [49]. Our results suggest that genes in this model were initially involved in development, while the oldest group exhibited a reduction of brain function; in other words, chicken cerebrum senescence was accompanied by impairment of neural plasticity, increases in inflammatory response and neurological disorders, and decrease of learning and memory abilities. Generally, hens at the age of 300 days are high egg productive, subsequently, their physiological features change substantially with age, leading to a reproductive decline and a drop in egg laying. Therefore, we can detect the biological variations in RNA levels.

Different complex temporal expression patterns

Genes in the green model entered periods of up-regulation during embryonic development (E12 to E16) and adulthood (D300 to Y1). On embryonic days 12 to 16, when histological differentiation is gradually completed, these up-regulated genes were enriched for the terms “myelin sheath” and “synaptic vesicle endocytosis”, indicating neural functions were important in this stage. Even so, the immune response was fully activated; for instance, the terms “lysosome”, “phagosome”, “complement activation”, “response to oxidative stress” and “T cell receptor signaling pathway” were enriched. It is well known that 300-day-old hens are adults, but we additionally detected signs of ageing at D300 to Y1, implying that the gene expression changes preceded physical ageing.

Functional prediction of lncRNAs

It is well known that lncRNAs are important regulatory factors exerted their roles by targeting corresponding mRNAs [50]. To gain insight into how putative lncRNAs performed the biological functions to regulate the cerebrum development and ageing, co-expression networks between DE lncRNAs and DEGs were constructed. Functional enrichment analysis of co-expressed genes showed similar results, except for pathways involved in the formation of neuron network, cell division and development, genes such as SYT1, CDH1, PAX6, NEFL, and ATXN, which are typical neuronal markers in brain development [2] were discovered in this study. Besides, GO terms and KEGG pathways related to defense mechanisms and immune responses were enriched, implying the biological function of DE lncRNAs in chicken cerebrum in regulating the brain developmental and immunological mechanisms.

Conclusions

In summary, we have identified a set of age-dependent lncRNAs and mRNAs, determined their temporal expression patterns and clarified the dynamic changes in the chicken cerebrum transcriptome over the course of time. Although some important genes and pathways were reported in this study, no work has been conducted on the functions of lncRNAs in chicken cerebrum, the molecular genetic mechanisms regulating chicken brain development and ageing are still poorly known. Therefore, to validate studies based on these hypothetical results deduced by bioinformatics methods were needed. Collectively, these findings may contribute to further studies on the molecular mechanisms of chicken development. In addition, as the chicken is an important model organism bridging the evolutionary gap between mammals and other vertebrates [51], our results may help to improve our comparative understanding of the cerebrum during vertebrate evolution.

Materials and methods

Ethics statement

The experimental procedures in this study were approved by the Committee on the Care and Use of Laboratory Animals of the State-level Animal Experimental Teaching Demonstration Center of Sichuan Agricultural University (Approval ID: B20160403). Anesthetization was performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (China, revised 2004), phenobarbital sodium solution with the concentration of 30 mg/kg was intraperitoneally injected and all efforts were made to minimize the suffering of the chickens. Finally, chickens were euthanized with overdose of phenobarbital sodium.

Chick embryo incubation and sample collection

Fertile eggs (Tibetan chicken) were bought from Mao County Jiuding Original Ecological Livestock and Poultry Breeding Co., Ltd., Aba Autonomous Prefecture, Sichuan Province. Fertile eggs were incubated at 37.8 °C (incubation days of 1–19) and 37.2 °C (incubation days of 20–21) with 70–80% humidity; a 24 h lighting schedule and intermittent rotations were employed. Ten chick embryos were randomly obtained on embryonic days 12, 16, and 20 (E12, E16 and E20). We sampled cerebrums for further high-throughput sequencing and skeletal muscle for sex determination. All chicks were housed in the same room from hatching to 6 weeks of age at a density of 12 chicks /m2 under uniform environmental conditions at Ya’an Poultry Farm, which belongs to Sichuan Agricultural University. After that, chicks were housed in cages, the uniformed poultry integration cage systems were established with 2 storeys, the size of each unit was set as 2.5 m length × 1.8 m width × 0.8 m height, and the density of stock was controlled to be 10 chicks /m2. All food, water, temperature, and light condition were monitored digitally. Healthy female hens were sacrificed at the age of 100 and 300 days (D100 and D300), together with hens aged 1, 3 and 5 years, cerebrum tissues were collected from 3 biological replicates at each time point except for E20, Y3 and Y5 (2 replicates). It was noted that, we collected 21 samples across 8 developmental stages in total, only 2 female embryos were identified among samples collected at the embryonic day 20, in order to collect samples in aging (Y5) as many as possible, the left hens were divided into 2 groups at Y3 and Y5 with 2 replicates, respectively. All the samples were snap frozen in liquid nitrogen immediately after collection and then stored at − 80 °C until RNA extraction.

Embryonic sex determination by PCR

The embryos were sexed by PCR, which was performed to amplify the CDH1-Z/W sequence. We designed CHD1-based primers because some introns of CHD1 differ between the Z and W chromosomes [52]. The sequence of the forward primer was 5′-GTTACTGATTCGTCTACGAGA-3′ and that of the reverse primer was 5′-ATTGAAATGATCCAGTGCTTG-3′. The PCR program was as follows: 5 min at 94 °C; 35 s at 94 °C, 30 s at 54 °C, and 30 s at 72 °C for 40 cycles; 5 min at 72 °C for final extension. Electrophoresis was performed on 1.5% agarose gel to identify the PCR products; males were identified by the presence of one band (CHD1-Z band), while females were identified by the presence of 2 bands (CHD1-Z and CHD1-W bands); female chicks were selected from every stage to use for high-throughput sequencing.

RNA isolation, library construction and Illumina deep sequencing

We used TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) to isolate total RNA following the manufacturer’s instructions; genomic DNA was removed by DNaseI (Qiagen, Beijing, China). The RNA purity was determined by measuring absorbance at 260 nm and 280 nm on an ND-1000 spectrophotometer (NanoDrop 2000, Thermo Fisher Scientific, Waltham, MA, USA) to calculate the A260/280 ratio; the integrity and concentration were estimated using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, Calif.) with the Agilent RNA Nano Kit. In addition, RNA integrity was examined by 1% agarose gel electrophoresis. Samples with RIN (RNA integrity number) values above 7.5 and 28S/18S ≥ 1.0 were used for library construction and sequencing. The Ribo-Zero™ Gold Kit (Illumina, San Diego, CA, USA) was used to remove rRNA from the total RNA. Subsequently, the sequencing libraries were generated following manufacturer recommendations with varied index label by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, USA). The details of library construction showed as follow: Firstly, rRNA was removed by kits, RNA fragmentation and short RNA strands were carried out by NEBNext First Strand Synthesis Reaction Buffer under elevated temperature. Subsequently,first cDNA strand was synthesized using random hexamer primers and RNA fragments as template. Second strand cDNA synthesis was subsequently performed using buffer, dNTPs, DNA polymerase I and RNase H. The library fragments were purified with QiaQuick PCR kits and elution with EB buffer, then terminal repair, poly (A)-tailing and adapter ligation were implemented. In order to select cDNA fragments of preferentially 300 bp in length, the library fragments were purified and the UNG enzyme was used to digest second strand of cDNA. PCR was performed by aiming expected size of amplicons, and the library was completed. The resulting 21 libraries were sequenced using the Illumina HiSeq platform with a paired-end sequencing length of 150 bp (PE150) at Annoroad Gene Technology Corporation (Beijing, PR China).

The clean data have been submitted to the NCBI Gene Expression Omnibus (GEO) with the accession number GSE114129.

Read mapping, quantification and transcriptome annotation

The quality control of raw data was completed by removing low-quality reads, adaptor sequences, empty reads and rRNA reads. Processed reads from each sample were mapped to the chicken reference genome (Gallus_gallus-5.0) using HISAT2 2.1.0 followed the previous described protocol [28, 53]. The mapped reads were assembled with StringTie v1.3.3 [54]; assembled transcripts were then merged into consensus transcripts using custom Python scripts [55] after filtering reads with length less than 200 nt. Finally, transcripts annotated as “c” and “=” were removed by Cuffcompare v2.2.1, noted that, “c” for partial match, “=” for full match [56]. Then, we obtained a set of putative coding and non-coding transcripts, and the known transcripts were removed by blastx and Hmmscan [57].

Identification and classification of lncRNAs

The Coding Potential Calculator (CPC) was applied to estimate the coding ability of the remaining transcripts, and transcripts with CPC score > 0 were filtered [58]. Subsequently, transcripts with FPKM > 0 at least in one biological replicate were identified as lncRNAs (Noted that FPKM is abbreviation of fragments per kilobase of transcript per million mapped reads). LncRNAs were initially classified into 2 major categories by FEELnc [59] according to their location with respect to mRNA: (1) intergenic lncRNAs, which do not intersect with any gene annotations; (2) intragenic lncRNAs, which overlap with gene annotations and can be further classified into 4 subcategories: (a) lncRNAs that overlap with sense introns, (b) lncRNAs that overlap with antisense introns, (c) lncRNAs that overlap with sense exons, and (d) lncRNAs that overlap with antisense exons. The intragenic lncRNAs occurred in two situations: either the lncRNA contained the intron or exon, or the lncRNA was contained in the intron or exon.

Transcriptomic gene expression analysis

The expression levels of the transcripts were expressed as FPKM values of mRNA and lncRNA using Cufflinks v2.2.1 [56] and StringTie v1.3.3 [54], where a value of 1 was added to the FPKM value of each gene before log2 transformation to avoid infinite values. FPKM > 0.1 was used to filter the expressed mRNA [60]. Pearson correlations were estimated across 8 developmental stages, and hierarchical clustering was performed using MultiExperimentViewer (MeV version 4.9.0) [61].

Age-dependent transcriptome and temporal expression pattern analysis

To test age-dependent differential gene expression at 8 time points, we identified differentially expressed mRNAs and lncRNAs using DESeq2 [62] based on the read counts produced from the FPKM calculated by StringTie v1.3.3 [54]. We conducted pairwise comparisons between adjacent time points. Genes and transcripts with FC ≥ 2 or ≤ 0.5 plus P-value ≤ 0.05 were identified as age-dependent genes. For instance, in the E16 vs. E12 comparison (the latter time point was defined as the numerator), a gene with the FC ≥ 2 would be grouped into the increased pattern and considered to be up-regulated during this stage, a gene with FC ≤ 0.5 would be grouped into the “decreased pattern”, and considered to be down-regulated. To evaluate the time course of age-dependent transcriptomic variations across the life cycle of the chicken, we performed STEM analysis (Short Time-series Expression Miner) based on differentially expressed genes and transcripts by comparing 2 adjacent time points [63]. Different colors were used to represent the significantly enriched model profiles (Bonferroni-adjusted P-values ≤ 0.05).

LncRNA-mRNA co-expression network analysis

We constructed a co-expression network using WGCNA (Weighted Gene Co-expression Network Analysis) in the R package [64] based on the differentially expressed mRNAs and lncRNAs that were identified by DESeq2 screening. Genes with coordinated expression patterns were termed “co-expression modules” and were detected by the dynamic tree cut method [65]. Functional annotation enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were conducted with the DAVID (Database for Annotation, Visualization, and Integrated Discovery) dataset (https://david.ncifcrf.gov) [66]. GO terms and KEGG pathways with P-values ≤ 0.05 were considered significantly enriched.

Validation of the gene expression by qPCR

We selected 10 DEGs and DE lncRNAs to validate the expression profiles of RNA-Seq by relative quantification/Real-Time PCR (q-PCR). The expression levels of the selected genes were normalized to β-actin, primers for the mRNAs/lncRNAs were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Additional file 1 Table S3). Total RNA was converted to cDNA using the EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen Biotech, China). The qPCR was performed on the Bio-Rad iQ5 Real-Time PCR Detection system to detect the RNA expression using TransStart Top Green qPCR SuperMix (Transgen Biotech, China). The reaction volume contained 5 μl TransStart Top Green qPCR SuperMix, 1 μl template cDNA, 0.4 μl forward and reverse primers, and added ddH2O to the final volume of 10 μl. Each sample was conducted with three technique replicates followed the qPCR system: 94 °C for 30 s, followed by 40 cycles of 94 °C for 5 s and 30s at the Tm indicated in Additional file 1 Table S3. The melting curve analysis was performed from 65 °C to 95 °C with increments of 0.5 °C for 5 s. The relative expression levels were calculated using the 2-△△Ct method, and the data were expressed as mean ± SD (standard deviation).

Abbreviations

DE lncRNAs:

Differentially expressed lncRNAs

DEGs:

Differentially expressed genes

FPKM:

Fragments per kilobase of transcript per million mapped reads

GEO:

Gene Expression Omnibus

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

STEM:

Short Time-series Expression Miner

WGCNA:

Weighted Gene Co-expression Network Analysis

References

  1. Carninci P, Kasukawa T, Katayama S, Gough J, Frith MC, Maeda N, Oyama R, Ravasi T, Lenhard B, Wells C. The transcriptional landscape of the mammalian genome. Science. 2005;309(5740):1559–63.

    Article  CAS  PubMed  Google Scholar 

  2. Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458(7235):223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, Devon K, Dewar K, Doyle M, Fitzhugh W. Initial sequencing and analysis of the human genome. Policies and Research in Identity Management. Springer Berlin Heidelberg, 2001.

  4. Birney E, Stamatoyannopoulos JA, Dutta A, Guigó R, Gingeras TR, Margulies EH, Weng Z, Snyder M, Dermitzakis ET, Thurman RE. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature. 2007;447(7146):799–816.

    Article  CAS  PubMed  Google Scholar 

  5. Sun M, Kraus WL. Minireview: long noncoding RNAs: new "links" between gene expression and cellular outcomes in endocrinology. Mol Endocrinol. 2013;27(9):1390–402.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Moran VA, Perera RJ, Khalil AM. Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res. 2012;40(14):6391–400.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Wang KC, Chang HY. Molecular mechanisms of long noncoding RNAs. Mol Cell. 2011;43(6):904–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Liu S, Wang Z, Chen D, Zhang B, Tian R, Wu J, Zhang Y, Xu K, Yang L, Cheng C. Annotation and cluster analysis of spatiotemporal- and sex-related lncRNA expression in rhesus macaque brain. Genome Res. 2017;27(9):1608.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, Fan L, Sandelin A, Rinn JL, Regev A. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res. 2012;22(3):577–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Telley L, Govindan S, Prados J, Stevant I, Nef S, Dermitzakis E, Dayer A, Jabaudon D. Sequential transcriptional waves direct the differentiation of newborn neurons in the mouse neocortex. Science. 2016;351(6280):1443.

    Article  CAS  PubMed  Google Scholar 

  11. Fatica A, Bozzoni I. Long non-coding RNAs: new players in cell differentiation and development. Nat Rev Genet. 2014;15(1):7–21.

    Article  CAS  PubMed  Google Scholar 

  12. Baumgart M, Groth M, Priebe S, Savino A, Testa G, Dix A, Ripa R, Spallotta F, Gaetano C, Ori M. RNA-seq of the aging brain in the short-lived fish N. furzeri - conserved pathways and novel genes associated with neurogenesis. Aging Cell. 2015;13(6):965–74.

    Article  CAS  Google Scholar 

  13. Lee CK, Weindruch R, Prolla TA. Gene-expression profile of the ageing brain in mice. Nat Genet. 2000;25(3):294–7.

    Article  CAS  PubMed  Google Scholar 

  14. Mercer TR, Dinger ME, Sunkin SM, Mehler MF, Mattick JS. Specific expression of long noncoding RNAs in the mouse brain. Proc Natl Acad Sci U S A. 2008;105(2):716–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Weindruch R, Prolla TA. Gene expression profile of the aging brain. Arch Neurol. 2002;59(11):1712–4.

    Article  PubMed  Google Scholar 

  16. Zeisel A, Muñozmanchado AB, Codeluppi S, Lönnerberg P, La MG, Juréus A, Marques S, Munguba H, He L, Betsholtz C. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science. 2015;347(6226):1138–42.

    Article  CAS  PubMed  Google Scholar 

  17. Bakken TE, Miller JA, Ding SL, Sunkin SM, Smith KA, Ng L, Szafer A, Dalley RA, Royall JJ, Lemon T. A comprehensive transcriptional map of primate brain development. Nature. 2016;535(7612):367.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Chodroff RA, Goodstadt L, Sirey TM, Oliver PL, Davies KE, Green ED, Molnã rZ, Ponting CP. Long noncoding RNA genes: conservation of sequence and brain expression among diverse amniotes. Genome Biol. 2010;11(7):1–16.

    Article  CAS  Google Scholar 

  19. Higaki S, Sato A, Kojima T, Oishi T. Microarray profiling of gene expression in the aging monkey brain. Neurosci Res. 2010;68(1):E305.

    Article  Google Scholar 

  20. Amaral PP, Mattick JS. Noncoding RNA in development. Mamm Genome. 2008;19(7–8):454–92.

    Article  CAS  PubMed  Google Scholar 

  21. Colantuoni C, Lipska BK, Ye T, Hyde TM, Tao R, Leek JT, Colantuoni EA, Elkahloun AG, Herman MM, Weinberger DR. Temporal dynamics and genetic control of transcription in the human prefrontal cortex. Nature. 2016;478(7370):519–23.

    Article  CAS  Google Scholar 

  22. Darmanis S, Sloan SA, Zhang Y, Enge M, Caneda C, Shuer LM, Hayden Gephart MG, Barres BA, Quake SR. A survey of human brain transcriptome diversity at the single cell level. Proc Natl Acad Sci U S A. 2015;112(23):7285–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Hawrylycz MJ, Beckmann CF, Ebbert A, et al. An anatomically comprehensive atlas of the adult human brain transcriptome. Nature. 2012;489(7416):391–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Kang HJ, Kawasawa YI, Cheng F, Zhu Y, Xu X, Li M, Sousa AM, Pletikos M, Meyer KA, Sedmak G. Spatio-temporal transcriptome of the human brain. Nature. 2011;478(7370):483.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Miller JA, Horvath S, Geschwind DH. Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci U S A. 2010;107(28):12698–703.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Rajagopal S, Pochet N, Srivastava G, Rosenkratz L, Mccabe C, Bennett D, Jager PD. Characterizing the role of MICRO-RNA expression in the aging human brain transcriptome. Alzheimers Dement. 2013;9(4):P828.

    Article  Google Scholar 

  27. Stern C. The chick: a great model system becomes even greater. Dev Cell. 2005;8(1):9–17.

    CAS  PubMed  Google Scholar 

  28. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Schug J, Schuller WP, Kappen C, Salbaum JM, Bucan M, Stoeckert CJ. Promoter features related to tissue specificity as measured by Shannon entropy. Genome Biol. 2005;6(4):R33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Wang Y, Xu HB, Wang MS, Otecko NO, Ye LQ, Wu DD, Zhang YP. Annotating long intergenic non-coding RNAs under artificial selection during chicken domestication. BMC Evol Biol. 2017;17(1):192.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Kuo RI, Tseng E, Eory L, Paton IR, Archibald AL, Burt DW. Normalized long read RNA sequencing in chicken reveals transcriptome complexity similar to human. BMC Genomics. 2017;18(1):323.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Hagelkruys A, Lagger S, Krahmer J, Leopoldi A, Artaker M, Pusch O, Zezula J, Weissmann S, Xie Y, Schöfer C. A single allele of Hdac2 but not Hdac1 is sufficient for normal mouse brain development in the absence of its paralog. Development. 2014;141(3):604–16.

    Article  CAS  PubMed  Google Scholar 

  33. Ma P, Schultz RM. HDAC1 and HDAC2 in mouse oocytes and preimplantation embryos: specificity versus compensation. Cell Death Differ. 2016;23(7):1119.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Ericson J, Rashbass P, Schedl A, Brennermorton S, Kawakami A, Van HV, Jessell TM, Briscoe J. Pax6 controls progenitor cell identity and neuronal fate in response to graded Shh signaling. Cell. 1997;90(1):169–80.

    Article  CAS  PubMed  Google Scholar 

  35. Marquardt T, Ashery-Padan R, Andrejewski N, Scardigli R, Guillemot F, Gruss P. Pax6 is required for the multipotent state of retinal progenitor cells. Cell. 2001;105(1):43.

    Article  CAS  PubMed  Google Scholar 

  36. Ruf F, Sealfon SC. Genomics view of gonadotrope signaling circuits. Trends Endocrinol Metab. 2004;15(7):331.

    Article  CAS  PubMed  Google Scholar 

  37. Naor Z, Benard O, Seger R. Activation of MAPK cascades by G-protein-coupled receptors: the case of gonadotropin-releasing hormone receptor. Cell. 2000;11(3):91–9.

    CAS  Google Scholar 

  38. Yang AL, Yang K, Meng-Ling WU: The relationship between polymorphism of FSHβ, POU1F1, FSHR gene and egg production traits in Tibetan chickens. J Southwest University Nationalities 2014.

  39. Chandrasekar B, Mummidi S, Perla RP, Bysani S, Dulin NO, Liu F, Melby PC. Fractalkine (CX3CL1) stimulated by nuclear factor kappaB (NF-kappaB)-dependent inflammatory signals induces aortic smooth muscle cell proliferation through an autocrine pathway. Biochem J. 2003;373(2):547–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Hughes PM, Botham MS, Frentzel S, Mir A, Perry VH. Expression of fractalkine (CX3CL1) and its receptor, CX3CR1, during acute and chronic inflammation in the rodent CNS. Glia. 2010;37(4):314–27.

    Article  Google Scholar 

  41. Sun H, Song J, Weng C, Xu J, Huang M, Huang Q, Sun R, Xiao W, Sun C. Association of decreased expression of the macrophage scavenger receptor MARCO with tumor progression and poor prognosis in human hepatocellular carcinoma. J Gastroenterol Hepatol. 2017;32(5):1107–14.

    Article  CAS  PubMed  Google Scholar 

  42. Mohamed A, Yang Z, Ning YY, Qin G, Raija S, Karl T, Lester K. The scavenger receptor MARCO is required for lung defense against pneumococcal pneumonia and inhaled particles. J Exp Med. 2004;200(2):267–72.

    Article  Google Scholar 

  43. Jonker MJ, Melis JP, Kuiper RV, Tv VDH, Wackers PF, Robinson J, Gt VDH, Dollé ME, Vijg J, Breit TM. Life spanning murine gene expression profiles in relation to chronological and pathological aging in multiple organs. Aging Cell. 2013;12(5):901.

    Article  CAS  PubMed  Google Scholar 

  44. Tozzini ET, Baumgart M, Battistoni G, Cellerino A. Adult neurogenesis in the short-lived teleost Nothobranchius furzeri: localization of neurogenic niches, molecular characterization and effects of aging. Aging Cell. 2012;11(2):241–51.

    Article  CAS  PubMed  Google Scholar 

  45. Somel M, Guo S, Fu N, Yan Z, Hu HY, Xu Y, Yuan Y, Ning Z, Hu Y, Menzel C, et al. MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain. Genome Res. 2010;20(9):1207–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Woll P, Morris J, Painschab M, Marcus R, Kohn A, Biechele T, Moon R, Kaufman D. Wnt signaling promotes hematoendothelial cell development from human embryonic stem cells. Blood. 2008;111(1):122.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Scheffzek K, Ahmadian MR. GTPase activating proteins: structural and functional insights 18 years after discovery. Cell Mol Life Sci. 2005;62(24):3014.

    Article  CAS  PubMed  Google Scholar 

  48. Chédotal A, Richards LJ. Wiring the brain: the biology of neuronal guidance. Cold Spring Harb Perspect Biol. 2010;2(6):a001917.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Gerrow K, Triller A. Synaptic stability and plasticity in a floating world. Curr Opin Neurobiol. 2010;20(5):631–9.

    Article  CAS  PubMed  Google Scholar 

  50. Pandey GK. Regulatory roles of noncoding RNA in development and disease; 2013.

    Google Scholar 

  51. International Chicken Genome Sequencing C. Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature. 2004;432(7018):695–716.

    Article  CAS  Google Scholar 

  52. Jarvi SI, Farias MEM. Molecular sexing and sources of CHD1-Z/W sequence variation in Hawaiian birds. Mol Ecol Notes. 2006;6(4):1003–5.

    Article  CAS  Google Scholar 

  53. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, Barrette TR, Prensner JR, Evans JR, Zhao S. The landscape of long noncoding RNAs in the human transcriptome. Nat Genet. 2015;47(3):199–208.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

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

  57. Eddy SR. A new generation of homology search tools based on probabilistic inference; 2009.

    Book  Google Scholar 

  58. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35:W345 Web Server issue.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, Jagannathan V, Cadieu E, David A, Lohi H. FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017;45(8):e57.

    CAS  PubMed  PubMed Central  Google Scholar 

  60. Consortium G. The genotype-tissue expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348(6235):648.

    Article  CAS  Google Scholar 

  61. Howe E, Holton K, Nair S, Schlauch D, Sinha R, Quackenbush J. MeV: MultiExperiment Viewer; 2010.

    Google Scholar 

  62. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. Ernst J, Barjoseph Z. STEM: a tool for the analysis of short time series gene expression data. Bmc Bioinformatics. 2006;7(1):191.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  64. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics. 2008;9(1):559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  65. Shi K, Bing ZT, Cao GQ, Guo L, Cao YN, Jiang HO, Zhang MX. Identify the signature genes for diagnose of uveal melanoma by weight gene co-expression network analysis. Int J Ophthalmol. 2015;8(2):269–74.

    PubMed  PubMed Central  Google Scholar 

  66. Aprea J, Prenninger S, Dori M, Ghosh T, Monasor LS, Wessendorf E, Zocher S, Massalini S, Alexopoulou D, Lesche M. Transcriptome sequencing during mouse brain development identifies long non-coding RNAs functionally involved in neurogenic commitment. EMBO J. 2013;32(24):3145–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We would like to thank the reviewers for their helpful comments on this article.

Funding

This study was supported by Fok Ying-Tong Education Foundation for Yong Teachers in the High Education Institutions of China (161026); the program from Sichuan Agricultural University (02920400), which were supported by the corresponding author Diyan Li.

Availability of data and materials

All the data supporting our findings are contained within the manuscript.

The clean data have been submitted to the NCBI Gene Expression Omnibus (GEO) with the accession number GSE114129.

Author information

Author notes

  1. Zhongxian Xu, Tiandong Che and Feng Li contributed equally to this work.

Authors and Affiliations

  1. Institute of Animal Genetics and Breeding, College of Animal Science and Technology, Sichuan Agricultural University, Chengdu, 611130, Sichuan, China

    Zhongxian Xu, Tiandong Che, Kai Tian, Qing Zhu, Shailendra Kumar Mishra, Mingzhou Li & Diyan Li

  2. Key Laboratory of Southwest China Wildlife Resources Conservation (Ministry of Education), China West Normal University, Nanchong, 637009, Sichuan, China

    Feng Li

  3. Novogene Bioinformatics Institute, Beijing, 100083, China

    Yifei Dai

Authors
  1. Zhongxian Xu
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tiandong Che
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Feng Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Kai Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Qing Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Shailendra Kumar Mishra
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Yifei Dai
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Mingzhou Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Diyan Li
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ZX performed the experiments, prepared figures and/or tables and wrote this manuscript. TC was responsible for all software used, analysis the data and prepared figures and Tables. FL performed the experiments, contributed materials and prepared the draft of this article for initial submission. KT performed the experiments, contributed materials and reviewed drafts of the paper. QZ contributed reagents/materials, reviewed drafts of the paper. SKM edited the manuscript and performed the verification experiments. YD performed the experiments and prepared the figures. ML conceived and designed the experiments, reviewed drafts of the paper. DL conceived and designed the experiments, and prepared the final response to reviewers and the revised article. All authors participated in writing and/or editing this manuscript and approved the final version of the manuscript.

Corresponding authors

Correspondence to Mingzhou Li or Diyan Li.

Ethics declarations

Ethics approval and consent to participate

The experimental procedures in this study were approved by the Committee on the Care and Use of Laboratory Animals of the State-level Animal Experimental Teaching Demonstration Center of Sichuan Agricultural University (Approval ID: B20160403).

Consent for publication

Not applicable.

Competing interests

The authors declare that 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. Data statistics of samples for RNA sequencing. Table S2. Number of mRNA and lncRNA (including exonic, intronic and intergenic lncRNA) transcripts distributed on chromosome identified in chicken cerebrum. Table S3. Primer sequences for qPCR. Figure S1. Pipeline of lncRNAs identification. Figure S2. Chromosome distribution of 3 subgroups of lncRNAs. The blue line represents the size of chromosome (the ratio of each chromosomal size to the total genomic size). The Pearson correlation of chromosomal size and exonic lncRNA, intronic lncRNA and intergenic lncRNA were 0.0960, 0.9789 and 0.9874, respectively. Figure S3. Age-dependent mRNAs (A) and lncRNAs (B). Age-dependent genes in each stage were identified through comparing 2 adjacent time points, genes and transcripts with the FC ≥ 2 or ≤ 0.5 plus the P-value ≤ 0.05 were identified as age-dependent genes. (DOC 295 kb)

Additional file 2:

Functional categories of genes overlapping with intronic and exonic lncRNAs. (XLSX 27 kb)

Additional file 3:

Top 10 significantly enriched categories of age-dependent genes in 7 clusters. (XLSX 21 kb)

Additional file 4:

Functional categories of genes based on significantly enriched models in STEM analysis. (XLSX 23 kb)

Additional file 5:

Functional categories of co-expressed genes in the top 2 modules. (XLSX 41 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

Xu, Z., Che, T., Li, F. et al. The temporal expression patterns of brain transcriptome during chicken development and ageing. BMC Genomics 19, 917 (2018). https://doi.org/10.1186/s12864-018-5301-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5301-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us