Biological functional annotation of retinoic acid alpha and beta in mouse liver based on genome-wide binding

Abstract

Retinoic acid (RA) has diverse biological effects. The liver stores vitamin A, generates RA, and expresses receptors for RA. The current study examines the hepatic binding profile of two RA receptor isoforms, RARA (RARα) and RARB (RARβ), in response to RA treatment in mouse livers. Our data uncovered 35,521, and 14,968 genomic bindings for RARA and RARB, respectively. Each expressed unique and common bindings, implying their redundant and specific roles. RARB has higher RA responsiveness than RARB. RA treatment generated 18,821 novel RARB bindings but only 14,798 of RARA bindings, compared with the control group. RAR frequently bound the consensus hormone response element [HRE; (A/G)G(G/T)TCA], which often contained the motifs assigned to SP1, GABPA, and FOXA2, suggesting potential interactions between those transcriptional factors. Functional annotation coupled with principle component analysis revealed that the function of RAR target genes were motif dependent. Taken together, the cistrome of RARA and RARB revealed their extensive biological roles in the mouse liver. RAR target genes are enriched in various biological processes. The hepatic RAR genome-wide binding data can help us understand the global molecular mechanisms underlying RAR and RA-mediated gene and pathway regulation.

Keywords: ChIP-Seq; cistrome; liver; motif; retinoic acid.

Fig. 1.

Hepatic genome-wide binding of retinoic acid receptors (RARA and RARB) in control and RA-treated mouse livers. A and B: total number of bindings and target genes, respectively, for RARA and RARB in response to RA. RARB has higher RA responsiveness than RARA. C and D: showed the number of overlapping and unique binding sites and target genes between 2 indicated groups. I: comparison between RARA and RARB in control. II: comparison between RARA and RARB after RA treatment. III: comparison between control and RA treatment in RARA bindings and target genes. IV: comparison between control and RA treatment in RARB bindings and target genes. The right and left percentage numbers displayed in the overlapped area represent the percentages of the overlapped peaks or target genes in right and left circles, respectively.

Fig. 2.

Analysis of DNA regions colocalized in RARA and RARB bindings in control and in response to RA. A: calculation of resolution (R), which is used to represent the relative location between 2 colocalized binding sites. Two bindings with R values of 0.5, 0.75, and 1 are shown as examples. B: distribution of resolutions for colocalized binding sites between the 2 indicated groups. “C” and “R” represent control and RA treatment while “a” and “b” represent RARA and RARB. C: distribution of fold enrichment ratios (log2) for colocalized bindings between the 2 indicated groups. D: number of binding sites per overlapped target gene. Each point represents 1 overlapped gene. Dashed lines represent the theoretical trend if the 2 indicated nuclear receptors have comparable binding density and solid lines represent the actual trend.

Fig. 3.

Comparison of hepatic bindings between RAR and RXRA in control and RA-treated mice. A: occupancy of RXRA in RAR bindings. B: number of genes uniquely targeted by each RAR or cotargeted by RAR-RXR. Percentages of RAR-RXR cotarget bindings and genes are shown in each column. C: distribution of fold enrichment ratios (log2) for colocalized bindings between RAR and RXRA (1-4). D: number of binding sites per overlapped target gene. Each point represents 1 overlapped gene between the 2 indicated nuclear receptors. Dashed lines represent the theoretical trend if the 2 indicated nuclear receptors have comparable binding density, and solid lines represent the actual trend. “C” and “R” represent “Control” and “RA treatment”; “a,” “b,” and “x” represent “RARA,” “RARB,” and “RXRA,” respectively.

Fig. 4.

Location of RARA and RARB binding sites in mouse genome. A: binding locations relative to the transcriptional start sites (TSS) of RAR target genes. B, C, and D: locations of RARA and RARB binding sites in control and RA-treated livers relative to the locations of the reference (ref.) genes included in Ref-Seq database. Chromosomes 1 (B), 12 (C), and X (D) were randomly selected as examples. Ca, Cb, Ra, and Rb represent the same definitions as indicated in the legend of Fig. 3. E: chromosomal distribution of RARB binding sites in RA-treated mouse liver. Each bar represents 1 binding site.

Fig. 5.

Motifs identified from RARA and RARB bindings in control and RA-treated mice. Based on P values given by the MEME-ChIP program, (A/G)G(G/T)TCA (A), CCC(G/A)CCCC (B), CCGG(A/G)A (C), and TAAAA(A/T) (D) are the 4 most common motifs in all RAR bindings (bottom). TOMTOM program identified the corresponding transcription factors to motifs shown in A, B, C, and D as nuclear hormone receptors, SP1, GABPA, and FOXA2, respectively (top). Heat maps showing co-occurrence of motifs for hormone response element (HRE), SP1 (SM), GABPA (GM), and FOXA2 (FM) in control RARA (E), RA-treated RARA (F), control RARB (G), and RA-treated RARB bindings (H). “SFM,” “SGM,” and “FGM” represent the binding with motifs of both SP1 and FOXA2, SP1, and GABPA, as well as FOXA2 and GABPA. “SFGM” represents bindings with 3 motifs (SP1, FOXA2, and GABPA) simultaneously. “Any” represents bindings with any 1 of the 4 motifs, while “None” represents bindings without any motifs. The number of bindings in each co-occurrence is color coded to the scale on the top right corner of each heat map.

Fig. 6.

Motif distribution and HRE repeat in RARA and RARB bindings. A1 and A2, B1 and B2, C1 and C2, and D1 and D2: distribution of HRE, SP1, GABPA, and FOXA2, respectively, in RARA and RARB bindings. E1–E4: representation of direct repeat profile in RAR bindings.

Fig. 7.

Profile differentiation for motif-associated RAR target genes at the biological process and pathway level by principle component analysis (PCA). Genes associated with a specific motif were grouped as a unit while the number of genes involved in various gene ontology biological processes (GOBP) or Kyoto Encyclopedia of Genes and Genome (KEGG) were treated as variables. A: PCA score plots based on data annotated with GOBP (left) and KEGG (right) database. Each shape and color represents a unit with a different combination of nuclear receptors, treatment groups (control as C and RA treatment as R), and motifs. Units with similar gene function profiles are grouped together within an ellipse. B: PCA loading plots based on data annotated with GOBP (left) and KEGG (right) database. Each point represents a biological process or KEGG pathway. C: identified GOBP (top) and KEGG pathways (bottom) that contribute significantly in differentiating the groups established in score plots. The numbers of genes involved in the various biological processes or pathways are color coded by the scale on the top right corner of each heat map. “C” and “R” represent control and RA treatment, while “a” and “b” represents RARA and RARB, respectively. “HRE,” “SP1M,” “GABPAM,” and “FOXA2M” represent HRE, SP1 motif, GABPA motif, and FOXA2 motif, respectively.

Fig. 8.

Identification of biological processes and pathways by orthogonal partial least square discriminant analysis (OPLSDA) to differentiate the biological function profile of 2 group of genes. Target genes of a specific nuclear receptor with the same motif were grouped as a unit while the number of genes involved in various biological processes or pathways were treated as variables. OPLSDA based on GOBP and KEGG database is shown in A1-A6 and B1-B6, respectively. A1-A6 and B1-B6, left: score plots. A1-A6 and B1-B6, right: S plots. The x-axis and y-axis of the score plots represent the scores of components 1 and 2. Based on the contribution (x-axis) and reliability (y-axis) shown in the S plots, up to 5 biological processes were identified to differentiate each 2 of the 4 groups shown in score plots of Fig. 7.

Fig. 9.

Visualization of the number of genes in biological processes and pathways identified from OPLSDA. The numbers of motif-associated genes involved in identified GOBP (A1-A6) and KEGG pathways (B1-B6) are color coded to the scale on the left of each heat map. “C” and “R” represent control and RA-treated mice, while “a” and “b” represent RARA and RARB.

Fig. 10.

Correlation between RA-dependent RAR bindings and hepatic gene expressions. A: number of RA-regulated hepatic genes that showed RA dependence in RARA and/or RARB binding. B: RA significantly induced the mRNA levels of Cyp26a1 and Rarb. C: RA induced RARB bindings in promoter region of Rarb gene (peaks 1 and 2). D: RA induced RARA and RARB bindings in the Cyp26a1 gene (peaks 1-8) *P < 0.01.