Overview of RNA transcriptomic profiles of chicken CAM tissues
After low-quality reads were removed, 110.3–114.1 M paired-end reads (2x150 bp in length) for the mRNAs and lncRNAs, and 23.8–29.2 M single-end reads (50 bp in length) for the miRNAs, were obtained from the CAM of Tibetan and Chahua chicken samples (Additional file 3: Table S3). The majority of the lncRNAs (60.7%) were intergenic long-stranded non-coding RNA (lincRNA), followed by intronic-lncRNA (17.7%), antisense-lncRNA (12.4%), and 9.2% sense-lncRNA (Fig. 1a). The length of mRNAs ranged from 51 to 30,0721 bp, while that of lncRNAs ranged from 202 to 223,626 bp; the majority of mRNAs and lncRNAs were 201–400 bp and >3000 bp long, respectively (Additional file 4: Fig. S1). In total, 20,912 expressed genes, 15,523 lncRNAs (3194 known), and 1439 miRNAs were identified in the CAM tissues of Tibetan and Chahua chickens. The 1439 miRNAs included 471 known and 968 novel miRNAs, and most of the miRNA reads were between 20 and 24 nt in length (Fig. 1b, Additional file 5: Table S4). The values of FPKM of mRNA/lncRNA and TPM of miRNA distribution were similar among the six samples (Fig. 1c, d, e).
Differential expression analysis of RNAs between TC and CH
We found 389 DE lncRNAs, 73 DE miRNAs, and 354 DE genes between the TC and CH groups, among which 194 DE lncRNAs, 46 DE miRNAs, and 180 DE genes were upregulated, whereas 195 DE lncRNAs, 27 DE miRNAs, and 174 DE genes were downregulated in the Tibetan chickens compared to Chahua chickens (Fig. 2, Additional file 6: Table S5). The DE lncRNAs, DE miRNAs, and DE genes were classified into two categories for the six samples, in which the expressed quantity of the DE lncRNAs, DE miRNAs, and DE genes had good repeatability within groups respectively (Additional file 4: Fig. S2).
Functional analysis of differentially expressed non-coding RNAs (DE lncRNAs and DE miRNAs)
For 389 DE lncRNAs, we predicted 722 cis- and trans-target genes that were mainly involved in the GO terms of angiogenesis, blood vessel development, regulation of vasoconstriction, regulation of blood pressure, response to oxygen levels, ATP metabolic process, fatty acid metabolic process, glucuronate metabolic process, enzyme binding, and ATP binding (Fig. 3a, Additional file 7: Table S6a). The KEGG analysis revealed the involvement of target genes in carbon metabolism, fatty acid metabolism, vascular smooth muscle contraction, and the Wnt, mTOR, MAPK, VEGF, and calcium signaling pathways (Fig. 3b, Additional file 7: Table S6b).
A total of 1766 target genes were predicted for 73 DE miRNAs. The target genes were mainly enriched in the GO terms of vasculogenesis, positive regulation of fatty acid oxidation, glycogen metabolic process, angiogenesis, vasculature development, response to hypoxia, oxygen metabolic process, vasodilation, ATP binding, and enzyme binding (Fig. 4a, Additional file 8: Table S7a). Notably, the enriched KEGG pathways were cysteine and methionine metabolism, amino sugar and nucleotide sugar metabolism, apoptosis, vascular smooth muscle contraction, fatty acid metabolism, carbon metabolism, vitamin B6 metabolism, and the MAPK, Notch, VEGF, mTOR, and calcium signaling pathways (Fig. 4b, Additional file 8: Table S7b).
Comparative parsing of lncRNAs and mRNAs and functional analysis of DE genes
To compare the exon number, open reading frame (ORF), and expression levels of lncRNA and mRNA, we needed to consider the differences in structure and sequence. The number of exons corresponding to lncRNAs was mainly less than 8, whereas the number of exons in mRNAs was relatively large, ranging from 1 to >30 (Additional file 4: Fig. S3a, b). The majority of lncRNA ORFs were mainly between 50 and 250 nt long, while the length of the ORFs of mRNAs was mainly between 100 and 1000 nt (Additional file 4: Fig. S3c, d). Expressed mRNAs and lncRNAs were located on all chromosomes of the chicken genome (Additional file 4: Fig. S4). The average and maximum values of mRNA expression were higher compared to those of lncRNA expression (Fig. 5a).
The 354 DE genes were mainly enriched in GO terms of defense response, homeostatic process, regulation of lipid metabolic process, glucose homeostasis, blood microparticle, oxygen transport, blood vessel endothelial cell migration, blood circulation, and respiratory system development. Enriched KEGG pathways included fatty acid metabolism, vascular smooth muscle contraction, as well as the calcium, PPAR and MAPK signaling pathways (Fig. 5b, Additional file 9: Table S8). Twelve DE genes (ACTC1, KCNMB4, NCS1, NGFR, ADAM8, CASQ2, IRF4, PTPRZ1, CALML3, ERBB4, ARR3, and NTSR1) were mainly associated with angiogenesis, blood circulation, hematopoiesis, as well as calcium and MAPK signaling pathways. Fifteen DE genes (SSTR5, NR1H4, HTR2C, APOA1, KCNB1, P3H2, CHST8, LYZ, HAO2, ACER1, ACSBG1, ELOVL2, ELOVL3, GBE, and NOX3) were related to glucose and carbohydrate metabolism, fatty acid metabolism, and lactate oxidation (Table 1).
Construction of the ceRNA network
We constructed a ceRNA network that yielded 529 pairs of candidate ceRNAs (lncRNA-miRNA-mRNA) from 162 DE lncRNAs and 108 DE genes through 25 miRNA target-mediated relationships (Additional file 4: Fig. S5). Considering that the network contains enormous information, and each relationship cannot be displayed in the figure, we constructed a mini-ceRNA (DE lncRNA-DE miRNA-DE mRNA) network of important RNAs. We identified 10 known mRNAs that were related to hypoxic adaptation and involved in angiogenesis and blood circulation (NGFR, ACTC1, CASQ2, ERBB4, KCNMB4, NCS1, and NTSR1) (Fig. 6a), as well as energy metabolism (SSTR5, NR1H4, and GBE) (Fig. 6b). A total of 39 lncRNA-miRNA-mRNA interactions were identified in the mini-ceRNA network constructed with these 10 DE genes and their targeted 37 DE lncRNAs and 9 DE miRNAs (Table 1, Additional file 10: Table S9).
Expression levels of DE lncRNAs, DE miRNAs, and DE genes using qRT-PCR
The expression levels of six DE lncRNAs (MSTRG.118362.24, MSTRG.25780.4, MSTRG.19949.13, MSTRG.29590.15, MSTRG.32082.2, and MSTRG.73129.28), six DE miRNAs (gga-miR-726-3p, gga-miR-6608-3p, gga-miR-6606-5p, novel-miR-676, gga-miR-460b-3p, and gga-miR-7442-5p), and seven DE genes (ARR3, NGFR, NTSR1, ADAM8, ACTC1, CASQ2, and CALML3) were measured using qRT-PCR to validate expression differences identified through RNA-seq. The results indicated that the expression levels of all RNAs were significantly different between TC and CH (P < 0.05), and fold-changes in expression followed the same trend in qRT-PCR and RNA-seq (Fig. 7).