Deep Sequencing Identification of Novel Glucocorticoid-Responsive miRNAs in Apoptotic Primary Lymphocytes

<div><p>Apoptosis of lymphocytes governs the response of the immune system to environmental stress and toxic insult. Signaling through the ubiquitously expressed glucocorticoid receptor, stress-induced glucocorticoid hormones induce apoptosis via mechanisms requiring altered gene expression. Several reports have detailed the changes in gene expression mediating glucocorticoid-induced apoptosis of lymphocytes. However, few studies have examined the role of non-coding miRNAs in this essential physiological process. Previously, using hybridization-based gene expression analysis and deep sequencing of small RNAs, we described the prevalent post-transcriptional repression of annotated miRNAs during glucocorticoid-induced apoptosis of lymphocytes. Here, we describe the development of a customized bioinformatics pipeline that facilitates the deep sequencing-mediated discovery of novel glucocorticoid-responsive miRNAs in apoptotic primary lymphocytes. This analysis identifies the potential presence of over 200 novel glucocorticoid-responsive miRNAs. We have validated the expression of two novel glucocorticoid-responsive miRNAs using small RNA-specific qPCR. Furthermore, through the use of Ingenuity Pathways Analysis (IPA) we determined that the putative targets of these novel validated miRNAs are predicted to regulate cell death processes. These findings identify two and predict the presence of additional novel glucocorticoid-responsive miRNAs in the rat transcriptome, suggesting a potential role for both annotated and novel miRNAs in glucocorticoid-induced apoptosis of lymphocytes.</p> </div

Pathways analysis predicts novel miRNA targets may contribute to glucocorticoid-induced apoptosis.

<div><p>(A) miRNA target predictions for novel miRNA candidates 44 and 166 were performed using the miRanda miRNA target prediction algorithm. The number of target mRNAs differentially expressed during glucocorticoid-induced apoptosis (p < 0.01; fold change > 1.2) is indicated for each candidate. </p> <p>(B) IPA-generated ranking of the top five molecular and cellular functions of genes differentially expressed during glucocorticoid-induced apoptosis (p < 0.01; fold change > 1.2), as well as the predicted targets of both candidates 44 and 166 (p-values for top functions are indicated beneath each ranking). Genes differentially expressed during glucocorticoid-induced apoptosis were identified by whole genome microarray analysis of untreated and 100nM dexamethasone-treated thymocytes (6 hours, 3 biological replicates). </p> <p>(C) Venn diagram analysis identified specific novel candidate predicted targets differentially expressed during glucocorticoid-induced apoptosis (p<.01) and the application IPA to this combined gene list (40 genes) generated a top 5 ranking of molecular and cellular functions regulated by these predicted targets (p-values for top functions are indicated beneath each ranking). </p></div

Development of a customized bioinformatics pipeline for the discovery of novel miRNAs from deep sequencing data.

<div><p>(A) This bioinformatics analysis workflow describes the novel miRNA discovery process adapted from miRanalyzer. The analysis pipeline uses next generation sequencing (miRNA-seq) data from untreated (control) or dexamethasone-treated rat primary thymocytes as input. This pipeline divides reads into three files: reads that align to an annotated mature miRNA (“Positive” training set), reads that align to other RNA subtypes (“Negative” training set), or reads that align at unannotated regions (“Test” set). Reads from each of these files are then aligned and alignment results are methodically processed to generate clusters, precursors and predicted secondary structures. Random forest machine learning is then employed to train the models for the prediction of novel miRNAs in the “Test” dataset. The output provides the genomic coordinates of predicted putative novel miRNAs.</p> <p>(B) Table describes total number of reads generated by miRNA-seq of control and dexamethasone treated primary thymocytes analyzed using the novel bioinformatics workflow described above. As expected, the majority of these reads align to known miRNAs when compared to other RNA subtypes. </p> <p>(C) Table summarizes the total number of known and predicted novel miRNAs identified by the bioinformatics workflow as induced or repressed in control and dexamethasone treated rat primary thymocytes. Both known and predicted novel miRNAs exhibit a trend of repressed expression during glucocorticoid-induced apoptosis.</p></div

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

<div><p>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.</p></div

E2f1 regulated and downstream pathways altered by AFB1.

<p>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.</p

Connection pathway analysis of DEGs from subchronic AFB1 exposure.

<p>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).</p

Novel exons found by RNA-Seq analysis.

<p>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.</p

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

<p>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.</p

Analysis for Differentially Expressed Genes (DEGs)^a.

a<p>Transcripts were grouped by various combinations of analysis into three columns: Total Assembled transcripts; RefSeq - transcripts matching or partially matching RefSeq genes; and Novel transcripts. In the first row, Total CuffCompare transcripts included all transcripts (57,076); RefSeq transcripts (complete or partial match) were 45,144 (14,257+30,887 = 45,144); and 11,932 potentially novel transcripts. Note, that a set of 1,496 Cufflinks assembled transcripts referred to as, ‘Others’, contained significant repeat sequences and so this small set of transcripts was excluded (e.g. –‘Others’) from Total Assembled transcripts. In row 2, to enable comparison with available Microarray data, we identified Cufflinks assembled transcripts that overlapped MA probes (total 44,469 transcripts). From this group of transcripts, we determined differential expression (DEGs) using student’s t-test in the case of MA (microarray) data in row 4, or for RNA-Seq data we performed DESeq analysis in row 3 or Cuffdiff analysis in row 5 (see <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0061768#s4" target="_blank">Materials and Methods</a> for details).</p

Methylation level of CpG island with an MBD3 peak that fall into within the indicated levels of DNA methylation as determined by RRBS.

<p>CpG islands bound by MBD3 (see <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004028#s4" target="_blank">Methods</a>) were binned by methylation level (see <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004028#s4" target="_blank">Methods</a>). Enrichment or depletion of MBD3 in each bin was determined by two-tailed t-test. Significantly enriched or depleted bins (p<0.001) are highlighted in bold.</p