RNA-Seq Profiling Reveals Novel Hepatic Gene Expression Pattern in Aflatoxin B1 Treated Rats

Merrick BA, Phadke DP, Auerbach SS, Mav D, Stiegelmeyer SM, Shah RR, Tice RR.
PLoS One (2013) DOI: https://doi.org/10.1371/journal.pone.0061768 PMID: 23630614

On This Page

Publication

Abstract

Deep sequencing was used to investigate the subchronic effects of 1 ppm aflatoxin B1 (AFB1), a potent hepatocarcinogen, on the male rat liver transcriptome prior to onset of histopathological lesions or tumors. We hypothesized RNA-Seq would reveal more differentially expressed genes (DEG) than microarray analysis, including low copy and novel transcripts related to AFB1's carcinogenic activity compared to feed controls (CTRL). Paired-end reads were mapped to the rat genome (Rn4) with TopHat and further analyzed by DESeq and Cufflinks-Cuffdiff pipelines to identify differentially expressed transcripts, new exons and unannotated transcripts. PCA and cluster analysis of DEGs showed clear separation between AFB1 and CTRL treatments and concordance among group replicates. qPCR of eight high and medium DEGs and three low DEGs showed good comparability among RNA-Seq and microarray transcripts. DESeq analysis identified 1,026 differentially expressed transcripts at greater than two-fold change (p<0.005) compared to 626 transcripts by microarray due to base pair resolution of transcripts by RNA-Seq, probe placement within transcripts or an absence of probes to detect novel transcripts, splice variants and exons. Pathway analysis among DEGs revealed signaling of Ahr, Nrf2, GSH, xenobiotic, cell cycle, extracellular matrix, and cell differentiation networks consistent with pathways leading to AFB1 carcinogenesis, including almost 200 upregulated transcripts controlled by E2f1-related pathways related to kinetochore structure, mitotic spindle assembly and tissue remodeling. We report 49 novel, differentially-expressed transcripts including confirmation by PCR-cloning of two unique, unannotated, hepatic AFB1-responsive transcripts (HAfT's) on chromosomes 1.q55 and 15.q11, overexpressed by 10 to 25-fold. Several potentially novel exons were found and exon refinements were made including AFB1 exon-specific induction of homologous family members, Ugt1a6 and Ugt1a7c. We find the rat transcriptome contains many previously unidentified, AFB1-responsive exons and transcripts supporting RNA-Seq's capabilities to provide new insights into AFB1-mediated gene expression leading to hepatocellular carcinoma.

Figures

Figure 1. DEGs identified from RNA-Seq by DESeq and Cuffdiff compared to microarray

Differentially expressed genes (DEGs) identified from RNA-Seq by DESeq and Cuffdiff compared to microarray.
Panel A. Number of DEGs observed with change ≥2-fold at p≤0.005 by DESeq, microarray and Cuffdiff analyses with Cufflinks assembled transcripts for all eight animals (4 Control, 4 AFB1) is shown. Total number of Cufflinks assembled transcripts and those that match an existing RefSeq annotation are displayed in the bar chart.
Panel B. Venn diagram of DEGs from Panel A shows the number of common transcripts (overlapping circles) and unique transcripts (non-overlapping circles) for all three analyses.
Panel C. Principal component (PC) analysis was performed for all samples using the gene expression values for DEGs found by microarray and DESeq (on Cufflinks assembled transcripts) analysis. The percentage variability captured by the first three principal components is displayed across PC#1, 2 and 3 represented on X, Y and Z axes.
Panel D. Heatmap shows all DEGs at ≥2-fold change, p≤0.005 from microarray or DESeq analyses (on Cufflinks assembled transcripts). Gene expression data were log2 transformed and then quantile normalized prior to generating the Heatmap for direct comparison of data. DEGs (red or green are upregulated or downregulated, respectively) for each animal were mapped by lane for each of four animals in the CNTL or AFB1 treatment groups

Figure 1 (6 MB)

Figure 2. Test and validation of differential expression between qPCR, RNA-Seq and microarray

Test and validation of differential expression between qPCR, RNA-Seq and microarray by select transcripts.
The mean fold changes for each group were compared in the bar chart for eight high to medium fold change transcripts and three low expression change transcripts (inset). qPCR data were normalized to β-actin expression for each sample. Means of significant (p≤0.05) fold changes from control were computed for qPCR, DESeq and microarray using RNA from the same 4 animals in each analysis

Figure 2 (804 KB)

Figure 3. Differing types of DEGs found by RNA-Seq and microarray data

Numbers of reads in RPM (reads per million mapped reads) are shown on the Y-axis and the genomic region is displayed on the X-axis using a common NCBI scale for representative AFB1 (blue reads) and CTRL (black reads) sample tracks in UCSC browser format. Placements of microarray probes for specific transcripts are indicated by small rose-colored squares with probe names below. Numbered Cufflinks assembled transcripts and the corresponding RefSeq (gene abbreviation) or Ensembl annotated transcripts (ENSRNOT) identifiers are displayed under the RNA-Seq tracks. Exons are represented as blocks or bands; introns are lines between exons. Arrows at the end of assembled or annotated transcripts show the direction of transcription. Open arrows in specific panels point out a new Cufflinks identified exon. Bar graphs to the right show mean fold changes (AFB1/CTRL) for specific microarray probes and the corresponding RNA-Seq transcripts (DESeq).
At the bottom of Panels A and C, the entire transcript (red) is presented for which a portion of the transcript which has been enlarged (blue bracket) for the RNASeq tracks shown above.
In Panel D, only spliced EST’s (parallel lines in the center of spliced EST represent gaps in the alignment) were available which were without RefSeq or Ensembl transcript annotation and for which no microarray probe had been assigned. In Panel D bar graph, fold changes for two separate Cufflinks transcripts are indicated while no microarray data were available (no probe).

Figure 3 (4 MB)

Figure 4. Novel transcripts found by RNA-Seq

Panel A. Bar graph shows the number of total and potentially novel transcripts on the Y-axis. Transcripts are shown for Total (purple), unique to either AFB1 (blue) or Control (red), and common (AFB1+ Control, green), according to the occurrence among replicates on the X-axis. A further breakdown of the 1811 total transcripts observed in two or more replicates is shown in the pie chart. Here, the majority of 1811 transcripts had some Ensembl annotation (1531) while the rest did not (280). Of those that were differentially expressed by AFB1, 21 had Ensembl annotation while 28 were unannotated and therefore, novel DEG transcripts.
Panels B and C show two such novel transcripts, described here as ‘HAfT-1’ and ‘HAfT-2’ (Hepatic Aflatoxin-responsive Transcripts 1 and 2). Numbers of reads as RPM are shown on the Y-axis and the genomic region is displayed on the X-axis for representative AFB1 (blue reads) and CTRL (black reads) sample tracks in UCSC browser format. See text for more details

Figure 4 (4 MB)

Figure 5. Novel exons found by RNA-Seq analysis

Panel A shows classification of different types of exons encountered during analysis of a Cufflinks assembled transcript in a model gene. Exons include ‘Exact’ matches to known exons, ‘Overlapping’ exons corresponding to partial matches, ‘Novel-T’ exons that occur with known transcripts, ‘Within’ exons occurring within the sequence of known exons and Novel-U exons that were unknown and occur outside known transcripts.
Panel B is a bar chart of such exon types that were unique to AFB1 or Control treatments and those exons that were shared between treatments.
Examples are shown for a Novel-T exon (Panel C) within the F11 transcript and two Novel-U exons (Panel D) outside the ASS1 gene. Numbers of reads, as RPM, are shown on the Y-axis and the genomic region is displayed on the X-axis for representative AFB1 (blue reads) and CTRL (black reads) sample tracks in UCSC browser format.

Figure 5 (4 MB)

Figure 6. Exon specific expression among homologous transcripts in the Ugt1a gene family

Panel A. The genomic region for Ugt1a transcripts is displayed on the X-axis in UCSC browser format where the Y-axis represents mapped reads in RPM units. Placements of microarray probes for specific transcripts are indicated by rose-colored boxes with probe names below. There are a total of four microarray probes, some of which correspond to shared exons of the Ugta1 gene family (A_44_P432355, A_44_P402641) or to specific exons defining Ugt1a1 (A_44_P446578) and Ugt1a6 (A_44_P432358) isoforms. RefSeq transcripts and Cufflinks assembled transcripts are displayed under the RNA-Seq tracks. Exons are shown as light blue blocks or bands; introns are lines between exons; arrows at the end of each transcript indicate direction of transcription.
Panel B. Bar graph shows mean fold changes (AFB1/Control) on the Y axis for the entire RNA-Seq transcript (blue), the microarray probe (red) and the isoform-specific RNA-Seq_exon (green). For some exons, there was no corresponding microarray probe − ‘No Probe’ (e.g. Ugt1a5, Ugta1a7c). Exon-specific, RNA-Seq ratios were labeled by exon number. Ugt1a-Common consists of four exons (common to all Ugt1a isoforms) for which two microarray probes exist. Exon-specific ratios from RNA-Seq reads were calculated for Exons 1, 2, 3 and 4. RNA-Seq exon-specific reads were measured to calculate AFB1/Control ratios for Ugt1a1, Utgt1a5, both exons of Ugt1a6, and Ugt1a7c.

Figure 6 (13 MB)

Figure 7. Connection pathway analysis of DEGs from subchronic AFB1 exposure

The top panel shows annotated interactions and regulatory relationships using IPA’s (Ingenuity Pathway Analysis) connectivity analysis. The connective pathway maps were generated using DEGs identified by DESeq_RNASeq (top panel) and for DEGs generated from microarray analysis (bottom panel) for only those transcripts with available RefSeq annotation. Hub genes (bolded, enlarged gene symbols) were defined as those transcripts regulating or interacting with ≥5 transcripts (red, upregulated; green, down-regulated).

Figure 7 (5 MB)

Figure 8. E2f1 regulated and downstream pathways altered by AFB1

AFB1 produced DEGs from DESeq analysis of RNA-Seq data which were analyzed by the IPA’s ‘grow pathway’ analysis (Ingenuity Pathway Analysis) which displays annotated regulatory relationships and interactions. Starting with induction of E2f1 in the center, DEGs from DESeq analysis were used to connect and grow downstream-dependent genes (red, upregulated; green, down-regulated). Hub genes (bolded, enlarged gene symbols) were defined as those transcripts regulating or interacting with ≥5 transcripts.

Figure 8 (4 MB)

Tables

Table 1. Alignment of RNA-Seq Reads to the Rat Genome

Table 1 (221 KB)