Integrated analysis of microRNAs, circular RNAs, long non-coding RNAs, and mRNAs revealed competing endogenous RNA networks involved in brown adipose tissue whitening in rabbits

Abstract

Background: The brown adipose tissue (BAT) is a target for treating obesity. BAT losses thermogenic capacity and gains a "white adipose tissue-like" phenotype ("BAT whitening") under thermoneutral environments, which is a potential factor causing a low curative effect in BAT-related obesity treatments. Circular RNAs (circRNAs) and long non-coding RNAs (lncRNAs) can act as competing endogenous RNAs (ceRNA) to mRNAs and function in various processes by sponging shared microRNAs (miRNAs). However, the roles of circRNA- and lncRNA-related ceRNA networks in regulating BAT whitening remain litter known.

Results: In this study, BATs were collected from rabbits at day0 (D0), D15, D85, and 2 years (Y2). MiRNA-seq was performed to investigate miRNA changes during BAT whitening. Then, a combined analysis of circRNA-seq and whole-transcriptome sequencing was used for circRNA assembly and quantification during BAT whitening. Our data showed that 1187 miRNAs and 6204 circRNAs were expressed in the samples, and many of which were identified as significantly changed during BAT whitening. Target prediction showed that D0-selective miRNAs were significantly enriched in the Ras, MAPK, and PI3K-Akt signaling pathways, and Y2-selective miRNAs were predicted to be involved in cell proliferation. The cyclization of several circRNAs could form novel response elements of key thermogenesis miRNAs at the back-splicing junction (BSJ) sites, and in combination with a dual-luciferase reporter assay confirmed the binding between the BSJ site of novel_circ_0013792 and ocu-miR-378-5p. CircRNAs and lncRNAs have high cooperativity in sponging miRNAs during BAT whitening. Both circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA triple networks were significantly involved in immune response-associated biological processes. The D15-selective circRNA might promote BAT whitening by increasing the expression of IDH2. The Y2-selective circRNA-related ceRNA network and lncRNA-related ceRNA network might regulate the formation of the WAT-like phenotype of BAT via MAPK and Ras signaling pathways, respectively.

Conclusions: Our work systematically revealed ceRNA networks during BAT whitening in rabbits and might provide new insight into BAT-based obesity treatments.

Keywords: BAT; Rabbits; Whitening; ceRNA; circRNA; miRNA.

Fig. 1

Identification and quantification of circRNA using circRNA-seq data and whole-transcriptome sequencing data. The transcripts marked by black lines and blue circles represent linear RNAs and circRNAs, respectively. The pooled RNA sample was used for the identification of circRNAs and the assembly of spliced circRNAs. The RNA samples in whole-transcriptome sequencing data were used to identify and quantify circRNA based on BSJ reads. The circRNAs identified by the two types of libraries were used for downstream analyses.

Fig. 2

Immunofluorescence (IF) assay and miRNA analysis of rabbit BATs at four growth stages. A IF of BATs at four stages. Two different magnifications are presented in the IF assay (200 μm and 20 μm). The nucleus were stained using DAPI. The blue and green signals represent the nucleus and UCP1 proteins, respectively. B Construction of circRNA- and lncRNA-related ceRNA networks during BAT whitening in rabbits. C Differential analysis of miRNAs. The scatter plot showed miRNAs with log2FC > 1 and < -1 in the six comparisons. The red and orange points show the miRNAs with a p-value < 0.05. D Validation of nine differentially expressed miRNAs (DEmiRNAs) using RT-qPCR. The red and black lines represent the miRNA-seq and RT-qPCR, respectively. The expression was normalized to the U6 gene and D0. The data shows the means of three independent experiments. Two technical replicates were set for one individual experimental replicate

Fig. 3

Target prediction and functional annotation of miRNAs. A Target genes of DEmiRNAs. The red and blue nodes show the miRNAs and mRNAs, respectively. Each edge shows the targeting relationship from miRNAs to mRNAs. B K-means clustering of DEmiRNAs based on expression levels. C Gene Ontology in biological process (GO-BP) enrichment and KEGG pathway analysis of target genes of miRNAs in miRC2. D GO-BP enrichment for the target genes of miRC4 and miRC7

Fig. 4

Genome-wide identification and differential analysis of circular RNAs (circRNAs) of BATs in rabbits. A Genomic distribution of circRNAs in a pooled RNA sample using circRNA-seq. B Venn diagram analysis of the circRNAs identified from circRNA-seq and whole-transcriptome data. C The Venn diagram analysis of expressed circRNAs among different growth stages. D Divergent primers amplify circRNAs. The IDs of circRNA were listed in Table S6. “M” represents “DNA marker”. The bands with expected length size were extracted and subjected to Sanger sequencing. The gels were cropped from Figure S5. E Sanger sequencing confirms head-to-tail splicing. The splicing mode and the GT/AG splicing signal of one circRNA were shown. The black arrows are pointing the back-splicing sites of the other eight circRNAs. F K-means clustering of DECs based on expression levels.

Fig. 5

Construction of circRNA-related ceRNA networks. A Construction of circRNA-related ceRNA network based on the head-to-tail linear sequences of circRNAs. The red, blue, and green nodes show the miRNAs, mRNAs, and circRNAs, respectively. B GO-BP enrichment and KEGG pathway analysis of circRNA-related ceRNA network. C An example of miRNAs interacting with the BSJ site of circRNAs. The black line shows the exons of MYCBP2. The 64th and 65th exon of MYCBP2 were highlighted and zoomed in a new window with a blue shadow. The predicted interaction between novel_circ_0013792 and ocu-miR-378-5p is shown in a yellow circle. The red block shows the BSJ site of circRNA. In the top-right panel, the interaction between miRNA and circRNA was shown in single base resolution. The black arrow points to the BSJ site of circRNA, which was located in the interaction region between miRNA and circRNA. The interaction between human hsa-circMYCBP2_107 and hsa-miR-378-5p is shown. D Prediction of the novel interaction pairs between circRNA BSJ sites and miRNAs

Fig. 6

Construction of lncRNA-related ceRNA networks. A Construction of lncRNA-related ceRNA network. The red, blue, and purple nodes show the miRNAs, mRNAs, and lncRNAs, respectively. B Comparing miRNAs involved in lncRNA-miRNA-mRNA networks and circRNA-miRNA-mRNA networks. C K-means clustering of lncRNAs based on expression levels. D GO-BP enrichment and KEGG pathway analysis of D85-selective lncRNA-related ceRNA network. E GO-BP enrichment and KEGG pathway analysis of Y2-selective lncRNA-related ceRNA network

Fig. 7

Validation of ocu-miR-378-5p-related ceRNA network. A Ocu-miR-378-5p-related ceRNA network. The green, purple, red, and blue nodes represent circRNA, lncRNA, miRNA, and mRNAs. B Dual-luciferase reporter assay validates the interaction between ocu-miR-378-5p and SLC15A3. The data in the dual-luciferase reporter assay shows the means of three independent experiments. The error bars show the standard errors of means. C Dual-luciferase reporter assay validates the interaction between ocu-miR-378-5p and MSTRG.16862.1. The data in the dual-luciferase reporter assay shows the means of three independent experiments. The error bars show the standard errors of means. D Validation of full length of novel_circ_0013792 by Sanger sequencing. “M” represents “DNA marker”. The bands with expected length size were extracted and subjected to Sanger sequencing. The gels were cropped according to the red lines in Figure S9. The black arrows point the BSJ sites of novel_circ_0013792 in Sanger sequencing. E Assembly of the full length of novel_circ_0013792. The blue and italicized bases show the sequence amplified by DIV1 primers and DIV2 primers, respectively. The black arrows show the binding sites and amplification directions of DIV1 primers. The red arrows show the binding sites and amplification directions of DIV2 primers (F) Dual-luciferase reporter assay validates the interaction between ocu-miR-378-5p and novel_circ_0013792. The data in the luciferase reporter assay shows the means of six independent experiments. The error bars show the standard errors of means