Genome-wide analysis of circular RNAs involved in Marek’s disease tumourigenesis in chickens

Figures & data

Figure 1. Identification of circRNA in chicken spleens. (a) The distribution of host genes encoding different number of novel circRNAs in chicken spleens. (b) The proportion of different categories of circRNAs in chicken spleens. Most of circRNAs were derived from exon regions of their host genes while a few circRNAs were formed by circularization of intronic or intergenic sequences. (c) The distribution of novel circRNAs in different chromosomes. Huge chromosomes produced more circRNAs in chicken genome. (d) Cumulative distribution of circRNA length (left) and expression (right). Three different types of circRNAs showed no significant changes in length and expression abundance (TPM). (e) Box plot showing the length of flanking introns of circRNAs’ host exon and the other exons. Introns surrounding circRNAs’ host exon significantly larger than that of other spliced exons. (f) Over-representation test showed that simple repeats and interspersed repeats were significantly enriched in circRNA flanking introns (p < 0.01, fisher exact test). Moreover, by employing alignment in each circRNA flanking intron pair, over 80% of RCMs contained chicken repeat 1 family (left). (g) The genomic regions of RUNX2. (h) The RCM in the long intron flanking sequences of the upstream of exon 2 and the downstream of exon 4.

Table 1. The specific information of 18 circRNAs.

Circ_ID	Chr	Parental gene	Length (nt)	Junction	Name
gga_circ_0011261	4	ZMYM3	272	Intron 10p^a	circZMYM3
ggs_circ_0014428	7	GTDC1	282	Intron 3p-exon 4^a	circGTDC1
gga_circ_0014634	7	MYO1B	1251	Exon 11–13	circMYO1B
gga_circ_0009435	3	RUNX2	546	Exon 2–4	circRUNX2.1^2
gga_circ_0009437	3	RUNX2	162	Exon 6	circRUNX2.2^b
gga_circ_0015132	8	AMY2A	605	Exon 5–9	circAMY2A
gga_circ_0003618	1	na	305	Intergenic_region	circNA1
gga_circ_0008071	2	UBE2E2	281	Exon 4–5	circUBE2E2
gga_circ_0002209	16	TAP1	315	Exon 6–7	circTAP1
gga_circ_0000062	10	na	181	Intergenic_region	circNA2
gga_circ_0000179	10	na	308	Intergenic_region	circNA3
gga_circ_0003057	1	RUNX1	411	Exon 2–3	circRUNX1
gga_circ_0000542	11	WWOX	357	Exon 3–4+ intron 3p^1	circWWOX
gga_circ_0003849	1	DACH1	278	Exon 2–3	circDACH1
gga_circ_0007434	2	na	507	Intergenic_region	circNA4
gga_circ_0016716	Z	NTRK2	718	Exon 2–3+ exon 10–11	circNTRK2
gga_circ_0016532	Z	SMARCA2	385	Exon 3–4	circSMARCA2
gga_circ_0005414	1	LRP6	329	Exon 8–9	circLRP6

^ap means the partial sequence of the region

^bcircRUNX2.1 and circRUNX2.2 were originated from the same parental gene but different exons.

Figure 2. Expressional profile of circRNAs. (a) Expression correlation between circRNAs and their respective parental genes. Reads that spanned junction site reads (JSRs) were compared between circular and linear RNAs. (b) Expression variance of circRNAs and their parental genes in three groups. (c) Heatmap showing the expression of all circRNAs identified in this research. Read counts of each circRNAs were scaled into Z score. (d) Number of differentially expressed circRNAs in three different contrasts. (e) Volcano plots depicting differentially expressed circRNAs in three contrasts. The horizontal dashed lines represent the adjusted p value cut-off of 0.01 and the perpendicular dashed lines represent |LFC| cut-off of 1.

Figure 3. GO enrichment of host genes of differentially expressed circRNAs. Directed acyclic graph (DAG) described the significant biological processes and molecular functions of DE circRNAs in TS vs. NS contrast (a) and TS vs. SS contrast (b).

Figure 4. Experimental validation of circular RNAs. (a) Divergent primers and convergent primers were used to amplify circRNAs in cDNA and genomic DNA (gDNA). Black triangles represent divergent primers white triangles represent convergent primers. (b) Sanger sequencing confirmed the back-splicing junction sequence of circRNAs.

Figure 5. Validation of the gene expression profile by real-time PCR. Six circRNAs were detected by Q-PCR. X-axis represents gene name, y-axis represents log₂foldchange, black column represents expression data from RNA-Seq, light grey column means expression data from Q-PCR.

Figure 6. The resistance of circRNAs to Rnase R digestion. (a) The expression of β-actin and parental genes of circRNAs with or without Rnase R treatment. (b) The expression of β-actin and circRNAs with or without Rnase R treatment. R+ represents that the RNA was digested by Rnase R, R- represents that RNA without digested by Rnase R. _**means p < 0.01.

Figure 7. Regulatory network of circRNA/lncRNA-miRNA-mRNA pairs in chicken spleens with Marek’s disease virus infection. (a) Interaction network of DE lncRNAs, circRNAs and mRNA. (b) Interaction network of miR-155 and its potential target lncRNAs, circRNAs and mRNAs.

Data availability

All data of this study are included in this published article and its supplementary information files. RNA-seq data are available at NCBI Gene Expression Omnibus (GEO, GSE138601 and GSE138600).

CircRNA data (GSE138600) are available at: (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138600), miRNA data (GSE138601) are available at: (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138601).