Differential Expression of Long Noncoding RNAs in the Livers of Female B₆C₃F₁ Mice Exposed to the Carcinogen Furan

Abstract

The mammalian genome is transcribed into mRNAs that code for protein and a broad spectrum of other noncoding (nc) RNA products. Long ncRNAs (lncRNA), defined as ncRNA species > 200 nucleotides long, are emerging as important epigenetic regulators of gene expression that are involved in a spectrum of biological processes of relevance to toxicology. We conducted a gene expression profiling study in the livers of female B₆C₃F₁ mice exposed to the carcinogen furan at 0.0, 1.0, and 2.0mg/kg (noncarcinogenic doses) and at 4.0 and 8.0mg/kg (carcinogenic doses) for 3 weeks. LncRNA differential expression showed a nonlinear dose response with none differentially expressed at 1.0 or 2.0mg/kg, 2 lncRNAs at 4.0mg/kg furan, and 83 at 8mg/kg, representing 13.3% (83/632) of the total number of differentially expressed transcripts. Among the lncRNAs observed, two lncRNAs examined showed transcriptional clustering with nearby protein-coding genes. LincRNA-p21 is an antisense transcript that is 15kb downstream from Cdkn1a locus and appears to be cotranscribed with the protein coding gene Cdkn1a at 8.0mg/kg furan. In a separate independent study, RNA samples from the livers of mice administered benzo(a)pyrene also demonstrated increased levels of Cdkn1a and the antisense lincRNA-p21 transcript. These data demonstrate that lncRNAs are transcriptional targets of furan exposures associated with levels of furan that are cytotoxic and induce cell proliferation. In addition, certain lncRNA transcripts are associated with the expression of nearby coding protein genes. We hypothesize that lncRNAs have potential as epigenetic biomarkers of carcinogenic exposures.

long noncoding RNAs, epigenetics, gene expression/regulation.

Although more than 60% of the content of the mammalian genome are transcribed, less than 2% of this is transcribed into mRNAs that code for proteins (Djebali et al., 2012). The transcriptional landscape of the mammalian genome was described in the study by Carninci et al. (2005) to consist of “a forest of transcription with overlapping transcription of the genome in both directions separated by deserts in which few transcripts are observed.” Interrogation of the RNA content of a variety of human tissues and cell lines has revealed that, collectively, at least 75% of genome is transcribed in one cell type or another (Clark et al., 2011; Djebali et al., 2012). Clark et al. (2011) analyzed the transcriptome of many mouse tissues and developmental stages and concluded that at least 63% of the genome is transcribed into thousands of protein-coding transcripts and over 30,000 long noncoding intronic, intergenic, and antisense transcripts. Although the pervasive transcription of the genome was once considered transcriptional noise, and the functions of most ncRNAs were unknown, it is apparent that certain ncRNA play critical roles as epigenetic regulators of gene expression for key biological processes of concern to the field of toxicology (Rinn and Chang, 2012).

Long ncRNAs (lncRNAs) are defined as RNA species > 200 nucleotides with no protein product. These RNA molecules were initially discovered by their association with transcriptionally active chromatin complexes (Guttman et al., 2009; Khalil et al., 2009; Rinn and Chang, 2012). The transcription of lncRNAs appears to be associated with their coorganization in the genome with transcriptional regulatory elements of protein coding genes. LncRNAs are frequently cotranscribed with nearby protein-coding genes from intergenic regions (referred to as long intergenic noncoding [linc] RNAs), whereas others are cotranscribed from within key protein-coding genes in the sense (lncRNAs) and antisense direction (referred to as natural antisense transcripts (NAT)) (Rinn and Chang, 2012)

To date, few lncRNAs have been functionally characterized (especially in rodents); however, those that have been characterized in humans are emerging as regulators of key cellular processes including stem cell pluripotency and differentiation (Guttman et al., 2009, 2011; Yang et al., 2011), cell cycle (Hung et al., 2011), molecular scaffolding (Guttman and Rinn, 2012), maintenance of DNA methylation/demethylation and transcriptional gene silencing (Mohammad et al., 2012), and global p53-mediated DNA damage response (Huarte et al., 2010; Yoon et al., 2012). Furthermore, certain lncRNAs are genetically imprinted and some lncRNAs are processed into functional microRNAs (Jalali et al., 2012; Moran et al., 2012). Functional studies have also demonstrated that lncRNAs fulfill critical roles as transcriptional and posttranscriptional regulators of gene expression by interacting with chromatin complexes (Guttman and Rinn, 2012; Nie et al., 2012). Once referred to as the “dark matter” or “junk” DNA of the genome, lncRNAs are emerging as novel epigenomic gene regulators that contribute to disease processes including cancer, neurodegenerative disorders, and diabetes (Gutschner and Diederichs, 2012; Wapinski and Chang, 2011). These data taken together suggest that lncRNAs can serve as integrated biomarkers of disease-related epigenetic effects.

In a furan dose-response study (0, 1.0, 2.0, 4.0, and 8.0mg/kg/day) using female B₆C₃F₁ mice (Moser et al., 2009), tumors were observed at 4.0 and 8.0mg/kg/day and were associated with hepatotoxicity and compensatory cell proliferation in a dose-related manner. Furan is considered to cause tumors in rodents by a nongenotoxic mode of action involving oxidative stress, accompanied by inflammation, cell proliferation, and toxicity (Ding et al., 2012; McDaniel et al., 2012). In this study, we conducted a short-term 3-week exposure of female mice to doses of 0.0, 1.0, 2.0, 4.0, and 8.0mg/kg/day furan and analyzed global expression profiles derived from liver tissue (Jackson et al., in preparation). This article is focused on the spectrum of lncRNAs that were induced by furan exposure in mice and the analysis by qRT-PCR of two specific lncRNAs: lincRNA-p21 and lncRNA Chromosome (Chr.) 9: 78107225–78118850 and nearby protein-coding genes. These two lincRNAs are coexpressed with nearby protein-coding mRNAs: lincRNA-p21 with Cdkn1a (cell cycle control gene, alias p21), and lncRNA Chr. 9: 78107225–78118850 with Gsta1 (encoding an enzyme that involved in detoxication), and Dppa5a (involved pluripotency/differentiation).

MATERIALS AND METHODS

Test article.

Furan (CAS No.110-00-9) used in this 3-week study (> 99% pure) was obtained from Sigma-Aldrich Chemical Co. (Milwaukee, WI). The appropriate quantity of the test article was mixed with Mazolas corn oil to a concentration that delivered the required dose levels. Each dose level was prepared separately on a volume-to-weight (vol:wt) ratio. Dosing solutions were dispensed into deaerated (with inert gas) 8-ml brown glass vials and capped and sealed with plastic closures adapted with silicon septa. Dosing solutions were stored in a refrigerator until use, but for no longer than 14 days. Previous studies have demonstrated the stability of furan dosing formulations for up to 14 days under these conditions (NTP, 1993).

Animals.

Female-specific pathogen-free (SPF) B₆C₃F₁ mice were received from Charles River Breeding Laboratories (Portage, ME) at 5–6 weeks of age and observed for approximately 7 days before study start. The animals were housed in cages in a SPF and AAALAC-accredited facility. All procedures were conducted in compliance with the Animal Welfare Act Regulations (9CFR1–4). Animals were handled and treated according to the guidelines provided in the National Institute of Health (NIH) Guide for the Care and Use of Laboratory Animals (ILAR, 1996). Feed (NIH-07; Zeigler Brothers, Inc., Gardners, PA) and tap water were available ad libitum. Mice were housed five per cage in polycarbonate cages.

Experimental design, necropsy, and liver sample collection.

Female B₆C₃F₁ mice were chosen for these studies because they were used in previous short-term studies with furan and were used in a furan cancer bioassay (Moser et al., 2009; Wilson et al., 1992). Male and female B₆C₃F₁ mice appear to respond similarly to furan-induced liver biochemical, proliferative, and cytotoxic effects (Fransson-Steen et al., 1997; Wilson et al., 1992). Female mice and the dose range of 1.0–8.0mg/kg body weight (bw) daily were chosen for this study because they represent no observed effect levels (1–2mg/kg) to dose levels (4.0–8.0mg/kg) known to induce a broad spectrum of biological effects including hepatocellular carcinomas at a dose of 8.0mg/kg.

Groups of five female B₆C₃F₁ mice received 0.0, 1.0, 2.0, 4.0, and 8.0mg furan/kg bw at a dosing volume of 10ml/kg bw by oral gavage daily for 3 weeks. On the final day of dosing, animals were euthanized 4h after final dose. Mice were anesthetized by CO₂ inhalation prior to euthanasia by exsanguination by severing the caudal vena cava after blood collection. The left, median, right posterior, and right anterior of the liver were cut into three sections. The right sections of each lobe were cut into 0.25–0.5cm³ pieces, combined from all lobes, frozen in liquid nitrogen, and placed into microcentrifuge tubes, snap frozen in liquid nitrogen and then stored at or below −70°C until RNA extraction.

RNA extraction from liver samples.

For RNA isolations, liver tissue was pulverized in liquid nitrogen and homogenized with a handheld Omni tissue homogenizer with an Omni generator probe (Omni International, Marietta, GA). Total RNA was extracted as per manufacturer’s instructions using the Qiagen RNeasy Midi Kit (Qiagen, Valencia, CA). The RNA was analyzed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA) for integrity, and concentration was determined using a Nanodrop ND-1000 (Thermo Fisher Scientific Inc., Wilmington, DE) UV-Vis spectrophotometer. Aliquots of total RNA were maintained at 4°C or stored at −70°C. Thawed aliquots of total RNA were shipped to Health Canada for microarray analysis.

Microarray analysis of furan liver samples.

The microarray analyses of protein-coding mRNAs and associated pathway and dose-response analyses will be described in a separate manuscript (Jackson et al., in preparation). RNA extracted from each mouse from all furan doses (0.0, 1.0, 2.0, 4.0, and 8.0mg/kg) were used for microarray analysis. Two hundred nanograms of sample RNA was used together with a mouse universal reference RNA (Stratagene by Agilent Technologies) to synthesize, amplify, and label cRNA using the Low Input Quick Amp Labeling Kit (Agilent Technologies). Labeled cRNA was purified using the RNease Mini Kit (Qiagen). Amplification and labeling efficiency of cRNA were quantified using a NanoDrop spectrophotometer. Hybridization mixes were prepared using the Hi-RPM Gene Expression Hybridization Kit (Agilent Technologies). Three hundred nanograms of Cy3-labeled reference RNA and 300ng Cy5-labeled sample cRNA were hybridized on SurePrint G3 Mouse GE 8x60K microarrays (Agilent Technologies) at 65ºC for 17h at 10rpm. Slides were washed according to the manufacturer’s specifications with Gene Expression Wash Buffers 1 and 2 (Agilent Technologies) and scanned using an Agilent G2505B scanner at 3 µm resolution. Data were extracted using Agilent Feature Extraction Software, version 11. The SurePrint G3 Mouse GE 8x60K microarrays are printed with 39,430 Entrez Gene RNAs and 16,251 lncRNAs.

Statistical analysis of microarray data.

Briefly, Agilent microarray data were preprocessed using Lowess with a log₂ transformation. Statistical analysis of gene expression changes was performed using ANOVA with contrasts between each dose and the vehicle control group. Probability values were adjusted for multiple comparisons using the “step up” False Discovery Rate (FDR) method of Benjamini and Hochberg (1995). Significant changes in gene expression were defined as p values < 0.05 and ± 1.5-fold change. The ANOVA analysis and associated contrasts were performed using Partek Genomics Suite (http://www.partek.com/). Pathway enrichment analysis was performed using both Partek and Ingenuity Systems (IPA) (http://www.ingenuity.com/). Significant enrichment p values were defined as a false discovery corrected p value < 0.05. The analysis of differentially expressed genes with respect to dose and pathways impacted by furan exposure will be published in a separate article (Jackson et al., in preparation). This article is focused on the finding that a high fraction of the probe sets differentially expressed in the furan-exposed animals was lincRNAs and the implication of this observation.

Quantitative RT-PCR.

Complementary DNA (cDNA) was generated from total RNA samples using the High Capacity cDNA Reverse Transcription Kit according to the manufacturer’s recommended protocol (Life Technologies, Carlsbad, CA). Quantitative real-time PCR (Q-PCR) was performed using TaqMan Fast Advanced Master Mix and custom TaqMan Gene Expression assays (Life Technologies) on the ViiA 7 Real-Time PCR System (Life Technologies): 50°C for 2min, 95°C for 20 s followed by 40 cycles of 95°C for 1 s and 60°C for 20 s. All Q-PCR primers used in the study were designed using Primer Express software (Life Technologies). The primer efficiency (E) of each custom TaqMan Gene Expression assay primer set was determined using the C_T Slope Method on serial dilutions of cDNA. Q-PCR analysis was performed using the ViiA 7 Real-Time PCR System software using the 2^−(∆∆C_T⁾ method.

Statistical analysis.

Data are presented as means ± standard error of the mean (SEM). For all endpoints, statistical analyses were done using Analyse-it statistical software (Analyse-it Software, Ltd, Leeds, UK). A one-way ANOVA was performed, and significant differences by ANOVA were further analyzed by Dunnett’s multiple comparison methods with p < 0.05.

qRT-PCR of lncRNA-p21 from livers of Benzo(a)pyrene-exposed mice.

The samples used in this analysis were from a previous experiment by Malik et al. (2012), and full experimental details are described therein. Briefly, adult Muta Mouse males were exposed for 28 days to 0, 25, 50, or 75mg/kg bw/day benzo(a)pyrene (BaP) dissolved in olive oil (n = 5 per group) by oral gavage. Three days after the last exposure, mice were euthanized by cardiac puncture under isofluorane anaesthesia. Mice were maintained under conditions approved by the Health Canada Animal Care Committee. Livers were removed, flash frozen, and stored at −80 until use. RNA was extracted and handled as described above. qRT-PCR was conducted and analyzed as described above.

RESULTS

Differential Expression of lincRNAs in Furan-Exposed Mice

Global expression profiling from the livers of female mice exposed to furan for 3 weeks showed that 2/11 (18.2%) at 4.0mg/kg and 83/632 (13.1%) at 8.0mg/kg of the differentially expressed probe sets on the SurePrint G3 Mouse GE 8x60K microarrays comprised lncRNAs (Table 1). No differentially expressed lincRNAs were observed at 1.0 or 2.0mg/kg furan in the liver. The lncRNAs in Table 1 are annotated by the Agilent probe set identification number and Chr. position number. The differentially expressed lincRNAs in mouse liver show lengths ranging from 199 to 80,751bp and are transcribed from the sense and antisense directions. With the exception of Chr. 16, each Chr. has 1–8 differentially expressed lncRNAs (Fig. 1). Only two lncRNAs among the 83 in Table 1 are annotated in the lncRNA database (http://www.lncrnadb.org/): lincRNA-p21 (Chr. 17: 29183003–29217681) and linc1257 (Chr. 11: 6421951–6429173). The chromosomal positions and lengths of other lncRNAs, as shown in Table 1, are from the Agilent database and the NCBI.

Table 1

Differentially Expressed lncRNAs at 4.0 and 8.0mg/kg Furan (p Value < 0.05 and ± 1.5-Fold Change*)

Agilent probeset ID	Chromosome number	lncRNA/chromosome genomic coordinates	lncRNA length (bp)	Strand	Fold change, 4.0mg/kg Furan	Fold change, 8.0mg/kg Furan
4.0 and 8.0mg/kg Furan
A_30_P01025511	4	lncRNA: 149445201–149445902	702	Forward	3.56	7.19
A_30_P01025159	6	lncRNA:129151512–129208737	57226	Reverse	−1.60	−2.07
8.0mg/kg Furan: 83 Unique lncRNAs
A_30_P01027205	1	lncRNA:138442536–138520766	78231	Forward	1.86	2.07
A_30_P01029956	1	lncRNA:163528200–163528398	199	Forward	1.66	1.89
A_30_P01026536	1	lncRNA:163508244–163586072	77829	Forward	1.30	1.64
A_30_P01026752	1	lncRNA:163528792–163529312	521	Forward	1.81	2.35
A_30_P01025359	1	lncRNA: 23340743–23341144	402	Forward	1.67	3.48
A_30_P01025912	1	lncRNA: 58443698–58449898	6201	Reverse	−1.43	−2.86
A_30_P01024134	2	lncRNA: 60038543–60046318	7776	Reverse	−1.16	−1.55
A_55_P2164075	2	lncRNA: 93292976–93310028	17053	Forward	−1.75	−5.72
A_30_P01026482	2	lncRNA:174945941–174946507	567	Reverse	−1.35	−1.55
A_30_P01018973	2	lncRNA:174943325–174966675	23351	Reverse	−1.37	−1.59
A_30_P01028205	2	lncRNA:175640275–175701875	61601	Forward	−1.30	−1.52
A_30_P01027279	2	lncRNA:176264400–176302925	38526	Forward	−1.26	−1.61
A_30_P01024689	3	lncRNA: 121615836–121616109	274	Forward	−1.62	−4.30
A_30_P01019901	3	lncRNA: 121627545–121639470	11926	Reverse	−1.05	−2.07
A_30_P01030305	3	lncRNA: 21964899–21974524	9626	Reverse	−1.41	−1.66
A_30_P01032887	3	lncRNA: 21973306–21974266	961	Reverse	−1.14	−1.58
A_30_P01027028	3	lncRNA: 51004595–51005002	408	Reverse	1.48	3.16
A_30_P01021547	4	lncRNA: 82156098–82196298	40201	Forward	−2.00	−–3.95
A_30_P01025368	4	lncRNA: 116791819–116798263	6445	Reverse	1.13	1.59
A_30_P01021813	4	lncRNA: 116790373–116797823	7451	Reverse	1.57	3.67
A_30_P01018287	4	lncRNA: 131734714–131741964	7251	Reverse	−1.10	1.56
A_30_P01023606	4	lncRNA: 146507925–146532650	24726	Reverse	1.23	1.68
A_30_P01025511	4	lncRNA: 149445201–149445902	702	Forward	3.56	7.19
A_52_P551011	4	lncRNA: 153394697–153400261	5565	Reverse	−1.48	−1.89
A_30_P01029305	5	lncRNA: 22887983–22939658	51676	Reverse	1.96	3.92
A_30_P01028374	5	lncRNA: 22890322–22939119	48798	Reverse	1.87	3.50
A_55_P2075258	5	lncRNA: 24607276–24615927	8652	Forward	−1.60	−2.18
A_55_P2361647	5	lncRNA: 25004885–25007359	2475	Forward	−1.35	−1.73
A_30_P01028938	5	lncRNA: 37051650–37110575	58926	Reverse	1.40	1.93
A_30_P01032429	5	lncRNA: 52461266–52541416	80151	Reverse	1.47	2.82
A_30_P01024631	5	lncRNA: 74487385–74497428	10044	Forward	1.25	1.53
A_30_P01025226	5	lncRNA: 121759042–121789992	30951	Reverse	1.61	2.29
A_51_P416613	6	lncRNA: 35202654–35213496	10843	Forward	−1.26	−1.71
A_30_P01017502	6	lncRNA: 108162092–108162886	795	Reverse	−1.27	−1.84
A_30_P01025159	6	lncRNA: 129151512–129208737	57226	Reverse	−1.60	−2.07
A_55_P2167955	7	lncRNA: 68427625–68456607	28983	Reverse	1.22	1.99
A_30_P01026270	7	lncRNA: 134620707–134626739	6033	Reverse	−1.97	−1.97
A_30_P01018981	7	lncRNA: 13533925–13544525	10601	Forward	2.58	3.77
A_30_P01030178	8	lncRNA: 27912000–27992750	80751	Forward	1.28	1.56
A_30_P01029190	8	lncRNA: 41475940–41476460	521	Forward	1.62	2.36
A_30_P01022473	8	lncRNA: 108120629–108123201	2573	Forward	−2.80	−2.21
A_30_P01029662	8	lncRNA: 10886200–10905925	19726	Forward	−1.70	−2.87
A_30_P01032012	8	lncRNA: 119739165–119739905	741	Forward	1.72	2.32
A_55_P2183513	8	lncRNA: 123113431–123117701	4271	Forward	−2.02	−3.36
A_30_P01031402	8	lncRNA: 124364048–124365705	1658	Forward	3.39	6.31
A_51_P185292	9	lncRNA: 34924313–34938507	14195	Forward	−1.54	−2.07
A_66_P104624	9	lncRNA: 72366069–72379665	13597	Reverse	1.81	2.47
A_30_P01021117	9	lncRNA: 75316023–75344401	28379	Forward	1.58	2.18
A_30_P01020518	9	lncRNA: 78104935–78116974	12040	Forward	1.40	4.00
A_30_P01018914	9	lncRNA: 78107225–78118850	11626	Forward	2.75	23.43
A_30_P01022249	10	lncRNA: 60473030–60535340	62311	Reverse	1.15	2.13
A_30_P01032909	11	lncRNA: 6421951–6429173	7223	Reverse	2.68	4.54
A_30_P01028002	11	lncRNA: 69189587–69190674	1088	Reverse	−2.40	−2.97
A_30_P01031674	11	lncRNA: 53596919–53672731	75813	Reverse	−1.69	−2.23
A_30_P01030240	11	lncRNA: 53592946–53637922	44977	Reverse	−1.45	−1.60
A_30_P01017718	11	lncRNA: 119797116–119804012	6897	Reverse	−1.22	−1.82
A_30_P01018231	12	lncRNA: 16833551–16901804	68254	Reverse	1.18	1.64
A_30_P01030946	13	lncRNA: 3497296–3500441	3146	Forward	1.54	2.62
A_30_P01024071	13	lncRNA: 65341237–65342006	770	Reverse	−2.59	−3.05
A_30_P01024541	13	lncRNA: 108259901–108272775	12875	Forward	1.34	2.15
A_30_P01018679	14	lncRNA: 73729977–73751877	21901	Reverse	1.65	2.50
A_30_P01031079	14	lncRNA: 73746324–73749119	2796	Reverse	1.72	2.28
A_30_P01020522	14	lncRNA: 65720392–65769859	49468	Reverse	−1.01	−1.96
A_30_P01021252	15	lncRNA: 77164450–77201325	36876	Forward	−1.26	−1.79
A_30_P01033325	15	lncRNA: 83266264–83294664	28401	Reverse	−1.50	−1.54
A_30_P01024684	15	lncRNA: 101084597–101095097	10501	Reverse	−2.19	−2.98
A_30_P01027464	17	lncRNA: 15100982–15115630	14649	Reverse	1.58	2.17
A_30_P01025984	17	lncRNA: 15226181–15250333	24153	Forward	2.24	2.41
A_30_P01021704	17	lncRNA: 17342991–17343357	367	Forward	−1.23	−1.53
A_30_P01020960	17	lncRNA: 29183003–29217681	34679	Reverse	2.21	5.43
A_30_P01023554	17	lncRNA: 34018371–34019394	1024	Reverse	1.24	1.73
A_55_P2075248	17	lncRNA: 35113844–35118733	4890	Reverse	−1.41	−1.64
A_30_P01032529	17	lncRNA: 85422370–85423168	799	Reverse	1.18	1.50
A_30_P01017543	17	lncRNA: 85418489–85427364	8876	Reverse	1.21	1.52
A_55_P2412319	18	lncRNA: 32518711–32537938	19228	Reverse	−1.25	−1.74
A_30_P01017688	18	lncRNA: 84303291–84303890	600	Forward	−1.37	−1.99
A_30_P01030779	18	lncRNA: 84757475–84757969	495	Reverse	−1.44	−1.60
A_30_P01024790	19	lncRNA: 23068397–23190256	121860	Reverse	1.13	1.99
A_30_P01018128	19	lncRNA: 36699712–36763966	64255	Reverse	−1.83	−2.80
A_30_P01026306	19	lncRNA: 57426515–57435165	8651	Reverse	−1.63	−2.07
A_30_P01021226	X	lncRNA:12811766–12813409	1644	Reverse	2.41	2.83
A_52_P223508	X	lncRNA: 47908921–47988498	79578	Reverse	−1.38	−1.64
A_55_P2310898	X	lncRNA: 148733043–148762038	28996	Reverse	−1.45	−2.49

Note. *Values shown in bold.

Fig. 1.

Chromosome distribution of differentially expressed lncRNAs based on microarray data in liver of furan-exposed mice.

LncRNAs and Clusters of Transcription

A number of lncRNAs show “clusters of transcription” with multiple transcripts originating from relative short segments of the genome: on mouse Chr. 1, three lncRNAs (199bp, 521bp, and 78kb) are transcribed from a 78kb chromosomal region; on mouse Chr. 2, four lncRNAs (567bp, 23kb, 62kb, and 39kb) are transcribed from a 1356kb chromosomal region; on mouse Chr. 4, four lncRNAs (702bp, 5.6kb, 7.3kb, and 25kb) are transcribed from a 21,665kb chromosomal region; and on Chr. 9, four lncRNAs (1.6kb, 12.0kb,13.6kb, and 28.3kb) are transcribed from a 5,752kb chromosomal region (Fig. 2A). Mouse Chr. 17 has three clusters of lncRNAs transcription (Fig. 2B). The “clusters of transcription” shown on Chr. 9 and 17 are nearby the protein-coding genes Gst1a (Chr. 9), Dppa5a (Chr. 9), and Cdkn1a (Chr. 17); these genes were also differentially expressed in furan-exposed mice (Jackson et al., in preparation). These three genes represent (1) Cdkn1a—a cell cycle gene known to be altered by a number of carcinogenic exposures in mouse liver (Harrill et al., 2009), (2) Gsta1—a member of the alpha class glutathione-S-transferases expressed in liver (Knight et al., 2007), and (3) Dppa5a (developmental pluripotency associated 5A)—a gene associated with stem cell pluripotency and differentiation (Kim et al., 2005).

Fig. 2.

Chromosome localization of differentially expressed lncRNAs and nearby protein-coding genes based on microarray data in liver of furan-exposed mice. (A) Chromosome 9; (B) Chromosome 17.

qRT-PCR of lncRNAs and Protein-Encoding Genes Cdkn1a, Gsta1, and Dppa5a

To validate the cotranscriptional regulation of specific clusters of lncRNAs/mRNAs by furan, we conducted qRT-PCR of a lncRNA (lncRNA Chr. 9:78107225–78118850) (Fig. 3A) nearby the protein-coding genes Gsta1 (Fig. 3B), Dppa5a (Fig. 3C), and lncRNA-p21 (Fig. 4A), an antisense lncRNA located 15kb downstream of Cdkn1a (Fig. 4B). We examined the degree of correlation of these RNAs in furan-treated livers. These three genes were among the top 10 genes impacted by furan exposure (Jackson et al., in preparation). qRT-PCR analysis revealed that the lncRNA on Chr. 9, in addition to the associated genes Gsta1 and Dppa5a, was not affected at 1, 2, or 4mg/kg furan relative to control mice. However, all three transcripts (lncRNA Chr. 9: 78107225–78118850, Gsta1, and Dppa5a) were significantly increased at 8mg/kg (p < 0.05) (Fig. 3). LincRNA-p21 and Cdkn1a mRNA each showed a dose response and were significantly elevated at 4 and 8mg/kg furan (p < 0.05) (Fig. 4). These data also suggest that exposure to furan induces bidirectional transcription at the protein-coding locus Cdkn1a and downstream in the antisense direction for lincRNA-p21 with both transcripts as potential independent targets for either p53-mediated transcription or a number of other transcription factors known to act at this site (Fig. 5) (Abbas and Dutta, 2009).

Fig. 3.

qRT-PCR of lincRNA Chr. 9 78107225–78118850 (A) and Gsta1 (B) and Dppa5a (C). Error bars indicate SE and *p < 0.05; 4.0 and 8.0mg/kg are carcinogenic exposure levels (Moser et al., 2009).

Fig. 4.

qRT-PCR of lincRNA-p21 (A) and Chr. 9 78107225–78118850 (A) and Cdkn1a (B). Error bars indicate SE and *p < 0.05; 4.0 and 8.0mg/kg are carcinogenic exposure levels (Moser et al., 2009).

Fig. 5.

The induction of a transcriptome response in furan-exposed mice induces the apparent bidirectional transcription of Cdkn1a and lincRNA-p21.

qRT-PCR of lncRNA-p21 in Livers of Mice Administered BaP

In order to determine whether the results obtained were specific to furan or reflected a broader response to hepatoxicants and hepatocarcinogens, we performed qRT-PCR for lincRNA-p21 on RNA isolated from the livers of mice exposed to the genotoxic carcinogen BaP from a previously reported study (Malik et al., 2012). Under the condition used in Malik et al. (2012), BaP induced a > 10-fold increase in Cdkn1a mRNA (Fig. 6A). RNA samples from this study analyzed by qRT-PCR also showed a > 5-fold induction of lncRNA-p21 (Fig. 6B), again suggesting bidirectional transcription of Cdkn1a and lncRNA-p21 in BaP-exposed mice.

Fig. 6.

qRT-PCR of Cdkn1a (A) (data reproduced from data in Malik et al., 2012) and lincRNA-p21 (B). Error bars indicate SE and *p < 0.05.

DISCUSSION

The demonstration that a major fraction of the eukaryotic genome is transcribed into tens of thousands lncRNAs has prompted the functional characterization of many of the ncRNA species. Although the functions of a large proportion of lncRNAs remain elusive, certain lncRNAs play key roles in biological processes and genome response networks that are known to be impacted by hepatotoxicants and hepatocarcinogens (Cui and Paules, 2010; Rinn and Chang, 2012). LncRNAs exert their biological effects through a number of mechanisms including direct binding to target mRNAs and interacting with chromatin regulators to control gene expression (Rinn and Chang, 2012). These current data suggest that lncRNAs are part of the cellular circuitry regulating differentiation, stress responses, and biochemical defense networks that can respond to cellular stressors (Rinn and Chang, 2012).

The observation that lncRNAs comprise a large part of the hepatocellular response to the liver carcinogen furan is novel. We hypothesized that it is not furan specific but rather reflects the transcriptome response to the biological impact of liver carcinogen-induced cytotoxicity and cellular proliferation. This global cellular defense response includes alterations in mRNA levels of protein-coding genes within cellular response pathways and as part of the epigenomic response network miRNA and lncRNAs. A recent article has shown the association with increased expression of an imprinted Dlk1-Dio3lncRNA cluster with hepatocyte hypertrophy in phenobarbital-exposed mice (Lempiäinen et al., 2013). The induction of the imprinted Dlk1-Dio3n lncRNA cluster was genetically dependent on both CAR and β-catenin signaling pathways, well-known pathways that impact mouse liver tumor development. A second expression profiling and mechanistic study using siRNA knockdown in livers of mice that had undergone 2/3 partial hepatectomy has also shown that a specific lncRNA (referred to as lncRNA associated with liver regeneration 1 or LALR1) has a key role in accelerating liver cell proliferation by activating the Wnt/β-Catenin signaling pathways (Xu et al., 2013). A previously published study from our group demonstrates dose-dependent induction of miR-34a and Cdkn1a mRNA in the liver of mice exposed to BaP (Malik et al., 2012). Reanalysis of RNA samples from this study by qRT-PCR also showed a dose-dependent induction of lincRNA-p21 (Fig. 6). These data taken together demonstrate that the transcriptional response to liver carcinogens includes not only protein-coding mRNAs and noncoding miRNAs but also includes lncRNAs. Moreover, the apparent pervasive transcription of the lncRNA epigenome in response to furan appears to be more robust than the limited miRNA epigenomic response of the liver to carcinogenic exposures such as BaP (Malik et al., 2012). These data demonstrate that the mode of action for the nongenotoxic carcinogen furan has an epigenetic component, and BaP has both a mutagenic and an epigenetic component to its mode of action.

Among the 83 lncRNAs showing differential expression in furan-exposed mice from this study, only 2 have been explored for functional properties in mouse and human cells: linc1257 and lincRNA-p21. Linc1257 is expressed in mouse embryonic stem cells and interacts with a number of chromatin-binding proteins/complexes in stem cells leading to repression of gene expression (Guttman et al., 2010, 2011). LincRNA-p21 is a repressor of p53-dependent transcriptional responses, affecting the expression of hundreds of p53 gene targets (Huarte et al., 2010). Based on studies using siRNA knockdowns, lncRNA-p21 appears to be necessary for p53-dependent induction of apoptosis but not cell cycle arrest (Huarte et al., 2010). This is mediated by physical association of lincRNA-p21 with heterogeneous nuclear ribonucleoprotein K (hnRNP-K), directing the genomic localization of hnRNP-K to repressed genes that trigger apoptosis (Yoon et al., 2012). These authors also demonstrated that lincRNA-p21 directly associates with Junb and Ctnnb1 mRNAs and selectively lowers their translation into proteins. These data indicate that lincRNA-p21 has two roles in regulating gene expression: (1) directly hybridizing to target mRNAs in the nucleus and (2) affecting the translation of p53 target genes in the cytoplasm.

Some of the lncRNAs observed in this study reside in the genome near protein-coding sequences and are apparently coexpressed with these genes. Furan resulted in the induction of a number of lncRNAs that demonstrated “clustering of transcription” with multiple lncRNA and protein-coding genes expressed within short chromosomal segments. We used qRT-PCR to validate two cases of apparent lncRNA coexpression with nearby protein-coding genes induced by furan exposure in liver. LincRNA-p21 is 15kb downstream from the start site of the p53 DNA damage-response gene Cdkn1a and appears to be coexpressed with Cdkn1a in the antisense direction (Fig. 5). In keeping with the response observed for furan, BaP exposure resulted in a highly similar dose-dependent increase in lincRNA-p21 (Fig. 6), mir-34a, and Cdkn1a expression in mouse liver (Malik et al., 2012). Studies using lincRNA-p21 siRNA knockdowns in hepatocytes cultures could be used to begin understanding impact of lincRNA-p21 on cancer hallmark biomarkers and the role of lincRNA-p21 in hepatotoxicity of carcinogenic exposures (Hanahan and Weinberg, 2011).

Numerous transcription factor binding sites are downstream of Cdkn1a and upstream of lncRNA including p53 (Abbas and Dutta, 2009), and in humans, this transcription factor–rich binding region is suggested to drive transcription of up to five lncRNAs from the Cdkn1a promoter in both the sense and antisense directions (Hung et al., 2011). In the human genome, two lncRNAs are downstream from the CDKN1 locus and are transcribed in the antisense direction: one of these is referred to as PANDA, whereas the second transcript is lincRNA-p21 (Hung et al., 2012). Both PANDA and lincRNA-p21 are part of the human cell p53 DNA damage response. In contrast, we have not shown increased levels for any other lncRNAs immediately downstream from Cdkn1a in the mouse. Bioinformatic examination of the downstream sequence of Cdkn1a also did not show a PANDA ortholog. The apparent species difference in the lincRNA-p21, PANDA, and Cdkn1a transcriptional cluster suggests a potential species difference in the p53 DNA damage response (Hung et al., 2012).

Recently, an experimental and computational survey of lncRNAs from three closely related rodent species, Mus musculus domesticus (Mmus: C57BL/6J), Mus musculus castaneus (Mcas), and Rattus norvegicus (Rnor), was conducted (Kutter et al., 2012). The authors used both RNA (RNAseq) and H3K4me3-bound (ChIPseq) DNA data (to identify putative active transcription sites) and combined both to construct catalogues of transcripts expressed in the adult liver of Mmus, Mcas, and Rnor. The authors found that a majority of the lncRNAs identified using this approach have no overlap with intergenic lncRNAs annotated in the mouse genome by Ensembl (build 64), demonstrating that current mouse lncRNA catalogues are largely incomplete. In addition, the lncRNAs observed showed mouse and rat strain differences in lncRNA transcripts; lineage-specific lncRNAs were associated with genomically neighboring protein-coding genes, and liver lncRNAs show differential expression when examined at five different developmental stages (e10, e12, e14, e18, and adult) of Mmus liver development. This study in mice and a previous study in humans have demonstrated the striking tissue specificity of lncRNA expression. For example, a survey of 4273 human lncRNAs revealed that 78% were tissue specific compared with only 19% of the 28,803 protein-coding genes (Cabili et al., 2011). Overall, these studies suggest that differential expression of lncRNAs may contribute to strain-, species-, and tissue-specific biological differences. Based on these findings, we also speculate that differential expression of lncRNAs may contribute to the spectrum of genomic and cellular responses to carcinogen exposure.

The purpose of this report was not to provide a comprehensive analysis of a toxicogenomic response of female mouse liver to furan, rather it was to demonstrate that lncRNAs are part of the transcriptional response to carcinogenic exposures, in this case furan. The robust increase of lncRNAs across the mouse genome in livers of female B₆C₃F₁ mice exposed to a cytotoxic carcinogen like furan is a novel observation, suggesting that lncRNAs may be epigenetic targets for mouse liver carcinogens. In addition, there is a clear need to assess if furan alters expression of lncRNAs in nontarget tissues.

The cellular processes influenced by lncRNAs can contribute to the six hallmarks of cancer that underlie malignant transformation (Gutschner and Diederichs, 2012; Hanahan and Weinberg, 2011). A survey of lncRNA expression within a panel of solid cancers has identified a number of novel transcribed regions differentially expressed across distinct cancer types that represent candidate biomarkers (Brunner et al., 2012). We propose that lncRNAs can serve as integrated epigenomic biomarkers of carcinogenic exposures since they can reflect a spectrum of biological processes of fundamental relevance to toxicology. Construction of an atlas of mRNA, miRNA, lncRNAs, and other mobile genetic elements from rodent and human tissues for use in toxicology studies is needed to assess the integrated transcriptional response of the genome to xenobiotic exposures as a tool to interrogate mode of action needed for biologically based risk assessments.

FUNDING

Work supported by ILS internal funding as indicated.

ACKNOWLEDGMENTS

The authors would like to thank the ILS Leadership of Dr TK Rao, Dr Thomas Goldsworthy, Dr David Allen, and Mr Sam Tetlow for funding this study and helpful discussions. We thank the ILS Genetic and Molecular Toxicology Division, Investigative Toxicology Division, and the Comparative Pathology Group for their contributions to this study. DNA microarray work was funded by Health Canada’s Genomics Research and Development Initiative. Salary for A.F.J. was provided through the Natural Sciences and Engineering Research Council of Canada.