Upregulation of the microRNA cluster at the Dlk1-Dio3 locus in lung adenocarcinoma.

Mice in which lung epithelial cells can be induced to express an oncogenic Kras(G12D) develop lung adenocarcinomas in a manner analogous to humans. A myriad of genetic changes accompany lung adenocarcinomas, many of which are poorly understood. To get a comprehensive understanding of both the transcriptional and post-transcriptional changes that accompany lung adenocarcinomas, we took an omics approach in profiling both the coding genes and the non-coding small RNAs in an induced mouse model of lung adenocarcinoma. RNAseq transcriptome analysis of Kras(G12D) tumors from F1 hybrid mice revealed features specific to tumor samples. This includes the repression of a network of GTPase-related genes (Prkg1, Gnao1 and Rgs9) in tumor samples and an enrichment of Apobec1-mediated cytosine to uridine RNA editing. Furthermore, analysis of known single-nucleotide polymorphisms revealed not only a change in expression of Cd22 but also that its expression became allele specific in tumors. The most salient finding, however, came from small RNA sequencing of the tumor samples, which revealed that a cluster of ∼53 microRNAs and mRNAs at the Dlk1-Dio3 locus on mouse chromosome 12qF1 was markedly and consistently increased in tumors. Activation of this locus occurred specifically in sorted tumor-originating cancer cells. Interestingly, the 12qF1 RNAs were repressed in cultured Kras(G12D) tumor cells but reactivated when transplanted in vivo. These microRNAs have been implicated in stem cell pleuripotency and proteins targeted by these microRNAs are involved in key pathways in cancer as well as embryogenesis. Taken together, our results strongly imply that these microRNAs represent key targets in unraveling the mechanism of lung oncogenesis.

Introduction

Lung cancer is the leading cause of cancer-related mortality in both men and women worldwide. Lung cancer subtypes include small cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC), which is subdivided into adenocarcinoma, squamous cell carcinoma and large cell carcinoma. Adenocarcinomas represent approximately 30% of all lung tumors. Each of these subtypes has different prognosis, disease signature and risk factors (e.g. smoking).

The role that KRAS plays during development, organ homeostasis and in the development of lung cancer has been extensively studied, leading to an understanding of its role in uncontrolled tumor growth, angiogenesis and inhibition of apoptosis.[1] KRAS is a GTP binding protein that is localized to the inner face of the plasma membrane. Guanine nucleotide exchange factors activate KRAS and enable cell growth and survival through downstream signaling pathways.[2, 3] KRAS is inactivated by GTPase activating proteins; however, many point mutations identified in KRAS prevent this GTP hydrolysis and thus maintain a constitutively active KRAS[4] KRAS mutations are very frequent in lung adenocarcinoma, occurring in approximately 25% of tumors[5, 6] Among KRAS mutations, variants at amino acid 12 represent 90 percent of the cases and a mouse bearing a conditional Kras mutation has been generated to study the effects of this mutation on cancer initiation.[5, 7, 8]

The decreased cost of sequencing technologies has sparked an interest in identifying the genetic changes that occur during tumor progression. Targeted re-sequencing[9] in addition to whole exome and whole transcriptome studies[10, 11] of tumor biopsies have provided several new candidate mutations. The analysis of the effect of additional somatic mutations and gene expression changes in mouse models of human tumors have complemented existing mutation data and provide a genetic framework for understanding tumor development. Additionally, genomic sequencing of several mouse strains has revealed several coding SNPs that can be used to identify parent of origin expression from F1 offspring from hybrid mouse strain crosses. Large scale non-coding RNA profiling studies have also identified microRNAs involved in oncogenesis including the miR-17 to miR-92 cluster of six microRNAs that are upregulated in B-cell lymphoma and SCLC[12, 13] and the common microRNA, let-7, has been shown to regulate the 3′UTR of Kras.[14] These previous studies have been insightful in understanding tumor development in lung adenocarcinoma, despite of focus solely on either coding or noncoding RNA populations.

We employed an integrative omics approach to identify transcriptional changes in a defined mouse model of lung adenocarcinoma. Sequencing was performed of bulk tumors, which reflect a mixed population of both tumor cells and infiltrating stromal cells; notable genetic changes could then be validated in purified sorted tumor cells. This approach revealed dysregulation of a complex stem cell associated microRNA locus in lung adenocarcinoma.

Results

Gene expression analysis indicates that a subset of genes are up and down regulated specifically in lung tumors

We performed a high throughput RNA sequencing analysis of the small and large RNA populations from three wildtype lungs and three Kras-driven lung adenocarcinomas. Two of these sets were derived from the offspring of an F1 between 129S4 and Molf/EiJ parents and one set from 129S4 homozygous parents (Table 1).

Table 1

Profile of samples subject to RNA sequencing.

	Tumor1	Tumor2	Tumor3	Normal1	Normal2	Normal3
Maternal strain	MOLF/EiJ	MOLF/EiJ	129S4	MOLF/EiJ	MOLF/EiJ	129S4
Paternal strain	129S4	129S4	129S4	129S4	129S4	129S4
Kras	G12D	G12D	G12D	wildtype	G12D*	wildtype
Cre-addition	Yes	Yes	Yes	Yes	No	No
p53 status	wildtype	wildtype	Inactive (flox/flox)	wildtype	wildtype	wildtype
Sequence reads	165,569,498	149,208,283	164,380,208	171,553,662	171,393,687	180,671,287
% mapped reads	87.17	84.62	86.01	80.76	84.14	81.45

Mutation is in inactive state until Cre recombinase is expressed

Differential expression analysis yielded ~450 significantly up-regulated and twice as many down-regulated genes in tumor versus normal lung samples (Figure 1a, Supplementary Table 1), distributed more or less evenly across the genome (Supplementary Figure 1). Genes significantly up in tumors include the Ros1 proto-oncogene (Supplementary Figure 2a) and Clec4n, both previously implicated in Kras transcriptomics and lung cancer.[15, 16] The top three genes decreased in tumors are the cGMP-dependent protein kinase Prkg1, the guanine nucleotide binding protein alpha (Gnao1) and the regulator of G-protein signaling 9 (Rgs9). These genes and Kras thus potentially are directly related or at least perform similar function as all are G-protein signal proteins. The Kras gene itself had similar levels of expression between tumors (p-value = 0.545; Supplementary Figure 2b).

Figure 1

Expression profile of aligned RNA sequences defines tumor versus normal differentially expressed sequences. (a) A scatterplot of mean FPKM values from all genes with a minimum FPKM of 0.01 (N = 21110). Genes significantly enriched in tumor samples are in green while those significantly down in tumor are red. Overall correlation between samples had an R2 value of 0.90782. (b) Gene expression profiling is comparable with microarray data from lung adenocarcinoma tumors.[17] Significant values from the microarray data set are plotted in the middle column. For genes significantly up in this data set, the corresponding RNAseq log2 fold change was plotted to the rightmost column. For genes significantly down in microarrays, the corresponding RNAseq values were plotted in the left-most column. (c) Venn Diagram of overlapping genes in this data set versus microarray examples.

Ingenuity Pathway Analysis of the most significantly differentially expressed genes provided broad categories of nucleic acid metabolism, embryonic/organ development and cell signaling and cancer that were enriched in differentially regulated genes in the tumor (Supplementary Figure 3a, 3b). In addition, pathway analysis revealed that direct and indirect connections could be established between the Kras gene, the p53 gene and the three most differentially down-regulated G-protein associated genes Prkg1, Rgs9 and Gnao1 (Supplementary Figure. 3c).

RNA expression profiles have been evaluated previously in similar Kras-driven lung tumors by microarray analysis that revealed 657 significantly differentially expressed genes.[17] The corresponding fold-changes of these genes in our RNAseq data were quite similar (Pearson r of 0.631; p < 0.001; Figure 1b) and several genes were considered significantly differentially expressed (112 up and 193 down) in each data set (Figure 1c).

Somatic mutations accumulate at similar frequencies in normal lung and tumor samples

The extensive coverage afforded by high throughput sequencing enabled us to identify 18, 34 and 23 nonsynonymous variants in the three tumor samples and 44, 40 and 80 variants in the three control lung samples (Supplementary Table 2). Genes with nonsynonymous variants from lung adenocarcinoma samples were not more frequently present in the Cosmic database[18] than the corresponding normal lung sample genes. Thus, in this model, Kras mutations do not appear to act together with multiple commonly mutated genes in lung or all cancers (Table 2).

Table 2

mRNA variant identification between samples.

	Tumor1	Tumor2	Tumor3	Normal1	Normal2	Normal3
GATK total variants	2,498,695	2,243,002	951,637	2,358,541	3,639,468	1,328,724
Passing filters	117,167	238,903	27,091	248,036	186,095	50,882
In exons	3,300	8,443	550	9,243	4,051	815
Unique	41	79	40	105	89	180
Synonymous	21	45	17	61	48	100
Nonynonymous	18	34	23	44	40	80
NS/S	0.86	0.76	1.35	0.72	0.83	0.80
Nonsyn genes in lung Cosmic mutations	5	0	1	6	5	9
Nonsyn genes in all Cosmic entries	11	28	10	42	32	73
Indels	13	23	6	19	23	14
Fusion transcripts	51	43	74	23	29	28

Analysis of known SNPs between MOLF and 129S4 mice identifies allele specific expression and potential areas of loss of heterozygosity

A total of 8065 coding variants in 4234 unique genes differentiate MOLF and 129S4 mice. This enabled analysis of allele specific expression of certain genes in addition to locations of loss of heterozygosity. Binning of SNPs based on their percent maternal expression did not show broad differences of parent-of-origin expression between tumor and normal lung samples (Figure 2a, Supplementary Figure 4). However, we could use the SNP information to identify individual genes with a biased allele specific expression (Supplementary Table 3). Cd22 had high levels in tumor samples (FPKM of 10.21 in tumor and 0.32 in normal lung; p=0.0001, Figure 2b). Surprisingly however, this expression largely or exclusively came from only the paternal allele while wildtype samples had bi-allelic expression (2 way ANOVA p-value of 0.028; Figure 2c). This primarily mono-allelic expression has been observed previously as a mechanism to retain specific antigen activity.[19] Sanger-based sequencing revealed no chromosomal amplifications and confirmed an allele-specific expression bias for Cd22 mRNA (Supplementary Figure 5). Allelic expression analysis for Kras revealed that neither the wildtype nor mutant allele was amplified, as can be the case in certain tumors with KRAS mutations.[20]

Figure 2

RNAseq analysis of coding variants reveals parent of origin specific expression or amplification of alleles. A total of 8065 coding variants that differ between 129S4/SvJae and MOLF/EiJ mice were interrogated for parent of origin expression. (a) Binning of percent maternal expression reveals that most SNPs follow a binomial distribution surrounding equal (50%) expression. (b) RNAseq read coverage for Cd22. Y-axis values represent read depth at each position adjusted by the total mapped reads for that sample relative to the mean mapped reads for all samples set from 0 to −300. Reads are negative because the gene is transcribed from right to left. (c) An analysis of genes that have the greatest fold difference of paternal to maternal allele expression reveals paternal-specific enhancement of Cd22 mRNA expression while maternal read counts remain similar between tumor and normal lung (2-way ANOVA p-value of 0.028 for tumor status). Error bars represent S.E.M. of the two F1 mice in each condition.

Analysis of RNA editing sites indicate that Apobec mediated C-to-U editing is common in tumors but not in wildtype lung

Post-transcriptional modifications, including RNA editing, can also be evaluated from RNAseq data. Several adenosine to inosine RNA editing sites have been identified in mice.[21] We found no evidence of differential adenosine to inosine editing in these tumor samples at these known positions (Figure 3a; Wilcoxon signed rank test p=0.19). However, the Apobec enzyme performs an alternate form of RNA editing, namely cytidine deamination leading to a uridine residue. Of the 30 editing sites that were identified in studies of Apobec1−/− mice[22] and were expressed in our lung samples, just over half (16) were C-to-U edited in tumor samples with editing ranging from 1.7% to 18.8% and the levels of editing were higher than corresponding rates in controls (Figure 3b). Detection of expression levels of Apobec1 revealed a modest 1.5 fold increase in expression in tumors (p<0.05; Figure 3c). A C-to-U edited site present in the Serinc1 3′UTR was validated by Sanger sequencing (Figure 3d). This data indicate that C-to-U editing is enriched in lung tumorigenesis, though it does not distinguish between whether the editing arises within the tumor during its formation or is the result of editing in immune derived cells in response to the tumor. Indeed, sorting for pure tumor originating cells (described below) revealed low Apobec1 mRNA levels and an absence of C-to-U editing at the Serinc1 site. While this cannot exclude the possibility that C-to-U editing in surrounding cells contributes to tumor progression it suggests that editing is not an inherent property of the tumor-originating cells.

Figure 3

Lung tumors display an increase in Apobec but not Adar RNA editing. (a) Mean percent editing for 168 previously identified Adar mediated A to I edited sites in the mouse with the size of the point reflecting the number of mapped reads at the given location. (b) Mean percent Apobec editing at 30 previously identified sites for the three samples, plotted as in panel A. (c) Apobec1 mRNA expression levels from the three sequenced mice in each condition plus tumors and normal lungs from four additional mice each, error bars represent SEM; * = p < 0.05 by t-test. (d) Sanger verification of C to T editing in the Serinc1 3′UTR. Primers were designed to amplify the edited site along with a coding SNP that indicated equal expression from both parental alleles.

Small RNA sequencing identifies a cluster of microRNAs upregulated in lung adenocarcinomas

Our most salient finding arose when we complemented RNA sequencing by performing small RNA sequencing on the normal lung and adenocarcinoma samples. Almost all of the most differentially expressed microRNAs aligned to an ~800kb region (nucleotides 110,691,433-111,519,307) on mouse chromosome 12qF1 (Figure 4a; Supplementary Table 4). 53 microRNAs align to this region with a mean fold increase of 31.6 (t-test p=3.8×10−87 relative to all other microRNAs; Figure 4b). This fraction of chr12qF1 microRNAs represents ~9% of all lung adenocarcinoma miRNA reads, compared with 0.1–0.2% of miRNA reads in the normal lung. Outside this locus the two closest microRNAs, miR-345 at ~600kb proximal and miR-203 at ~1.8mb distal, had equivalent expression between tumors and normal lung samples. This level of microRNA induction was substantially higher and more prominent than was observed for large RNAs suggesting that microRNAs may be key mediators of oncogenic drive in this mouse model.

Figure 4

A cluster of microRNAs on chromosome 12qF1 is upregulated in lung adenocarcinomas. (a) A scatterplot of microRNA counts normalized to one million mapped microRNA reads. Points in purple represent microRNAs that arise on chromosome 12qF1. (b) Boxplot of log2 based fold change for chr12qF1 microRNAs (in purple) and all remaining microRNAs (in black). (c) Small RNA northern blot of two representative microRNAs that align on chromosome 12qF1 (miR-127 and miR-376c), re-probed for tRNA sequences and/or U6. Small cell lung cancer (SCLC) and hepatocellular carcinoma (HCC) tumors were additional controls. (d) Cell lines derived from lung and liver tumor did not show high levels of chr12qF1 microRNA expression. (e) Quantitative RT-PCR analysis indicates that non-coding RNA expression of genes in the chromosome 12qF1 interval are highly expressed in lung adenocarcinomas relative to normal lung, SCLC tumors, and human lung adenocarcinoma tumor-derived cell lines. (* Probes for Rtl1 detect both sense and antisense transcripts) Values are plotted as mean +/− standard error of the mean of at least three samples run in triplicate. Note the differences in the y-axis for each of the plots. (f) microRNA qPCR indicates that an upregulation of the chr12qF1 locus is specific to tumors in vivo. Levels of two microRNAs from this cluster, miR-127-5p and miR-376a were normalized to U6 snRNA and then to wildtype lung levels. Cell lines derived from tumors were grown in culture or transplanted back into mice.

Northern blot analysis validated representative chr12qF1 microRNA expression patterns in lung adenocarcinoma samples from all 12 mice tested (two shown in Figure 4c). However, cell lines derived from the Kras mouse tumors, normal lung and SCLC samples had levels of chr12qF1 microRNAs comparable to wildtype lung (Figure 4c, 4d). Thus, the up-regulated microRNA cluster is a hallmark only of the tumors in vivo and not of the associated derived cell lines.

The mouse chr12qF1 region, also known as the Dlk1-Dio3 locus,[23] it is an area that is extensively methylated with a set of genes that are expressed specifically from the maternal or paternal chromosome. The cluster of microRNAs is transcribed from the chromosome inherited from the mother, as are the noncoding RNAs Meg3, Meg8, Rtl1as and Rian. Conversely, paternally expressed genes include Dlk1, Rtl1 and Dio3. SNPs in these genes were used to verify that this established parent of origin expression pattern was maintained in the tumor samples. From the RNAseq dataset, only the Rian gene was significantly up-regulated in the lung adenocarcinoma samples (FPKM of 6.48 for tumors versus 1.58; q-value < 0.005). However, when we examined previous microarray data from lung adenocarcinoma mouse tumors, probes for Meg3 were ranked second and sixth as the most abundantly represented in adenocarcinoma versus control (10.56 fold and 5.33 fold differences respectively) and a probe for Dlk1 ranked ninth on this list (5.18 fold increase).[17] We confirmed this up-regulation via quantitative real-time PCR with probes against Rtl1 (sense and antisense), Dlk1 and Meg3 (Figure 4e). Consistent with data from microRNAs in this locus, tumor samples displayed a marked increase in expression of these co-expressed noncoding RNAs (Figure 4e). Differentially methylated regions exist at the Meg3 promoter and intergenic to Meg3 and Dlk1 that are fully methylated in the adult lung. While hypomethylation at the maternal allele is typically associated with microRNA expression,[24] we observed no change in methylation in bulk tumors (Supplementary Figure 6) suggesting a different mode of locus activation is involved.

Increase of chr12qF1 microRNA expression is present in tumor-originating cells

Differential expression of genes or small RNAs from bulk tumors could arise from stromal cells that have infiltrated and become mixed with bona fide cancer cells. To reconcile whether the chromosome 12qF1 microRNAs and mRNAs are specific to the cancer cells, we generated tumors in Kras mice that also carried a Rosa26LSL-tdTomato Cre reporter allele. In this situation, when tumors were initiated with viral Cre, all resulting cancer cells would express red fluorescent protein along with the activated Kras allele. Fluorescence activated cell sorting of Tomato-positive cells from tumors enabled an evaluation of genes that are specific to the cancer cells. Validation of chromosome 12qF1 microRNAs revealed a ~35 fold induction of miR-127 and miR-376a in Kras expressing cells relative to wildtype lung (Figure 4f). Further, the Meg3, Dlk1 and Rtl1 genes were dramatically up-regulated (Supplementary Figure 7a). This demonstrates that the chromosome 12qF1 locus is activated in a manner specific to cells carrying a Kras mutation.

As noted above, cell lines derived from lung adenocarcinomas exhibited low levels of chr12qF1 gene and microRNA activation. However, when the Kras cells were transplanted back into mice, tumors arose within 3 weeks. By evaluating three of the tumors derived from transplanted cells, we were able to demonstrate that the same population of cells could reactivate their chr12qF1 microRNAs (Figure 4f) and mRNAs (Supplementary Figure 7b) when propagated in vivo. Taken together, this implies that the cell lines derived from lung adenocarcinomas do not display some of the more prominent features of solid tumors that are dependent on their in vivo environment.

miRNAs alter select target protein levels without globally influencing mRNA target expression

The large scale small RNA and large RNA sequencing datasets enabled global comparisons between predicted microRNA targets and their expression changes. The 53 microRNAs that are up-regulated on chr12qF1 were compared to 51 microRNAs that had a similar expression profile (within 10% of each other) between tumor and normal lung. Target mRNAs had no significant difference in expression for chr12qF1 microRNAs relative to the control microRNA set for Miranda (p=0.18), Pictar (p=0.46) and TargetScan (p=0.39) prediction programs (Figure 5a). This indicates that chr12qF1 microRNAs do not lead to reduction in mRNA levels of predicted targets.

Figure 5

Chr12qF1 microRNAs can repress protein levels of a subset of genes involved in oncogenesis. (a) The cumulative distribution of mRNA expression is unchanged for mRNAs predicted to be targets of the 53 chr12qF1 microRNAs upregulated in lung adenocarcinoma (in red) versus 51 control unchanged microRNAs (in blue) and overall mRNAs (in green). Three prediction programs were queried, Miranda (395 chr12qF1 targets and 376 control targets), Pictar (1231 and 1117 respective targets) and Targetscan (1002 and 1125 respective targets). (b) Luciferase expression of UTRs of mRNAs predicted to be targets of chr12qF1 microRNAs. Expression is normalized to firefly luciferase within the same construct and to a scrambled control shRNA (scr; black bars), transfected in E10.5 mouse embryonic fibroblasts (MEFs) where chr12qF1 microRNAs are not expressed. Significance was determined by a two-tailed t-test compared with a corresponding control shRNA (* = p < 0.05; ** = p < 0.01; *** = p < 0.001). Values represent the mean +/− SEM of at least two experiments performed in triplicate. (c) Site directed mutagenesis of three nucleotides of the miR-134 minding site in the Hyal1 3′UTR abrogated the repression of miR-134 on this 3′UTR. (d) Quantitative RT-PCR showed a reduction in Cftr but not Antxr1 mRNA levels in lung adenocarcinoma samples (n = 7) relative to normal lung (n = 7).

We utilized a luciferase reporter system to identify whether protein levels were altered by the up-regulation of these microRNAs. Twelve 3′UTRs that were predicted by TargetScan to be strong targets of ten of the most highly expressed chr12qF1 microRNAs. Co-transfection of the luciferase constructs and its corresponding target miRNA led to the down-regulation of 8 of these 11 targets from nine of the ten microRNAs relative to a scrambled control shRNA (Figure 5b). Three point mutations in the miR-134 binding site of the Hyal1 3′UTR caused a loss of miR-134 mediated repression of this 3′UTR (Figure 5c). Quantitative RT-PCR analysis revealed that Cftr but not Antxr1 showed a reduction in mRNA levels (Figure 5d). This is just a select subset of genes that are potential targets, yet it indicates that the microRNAs do have the intended effect of mediating post-transcriptional effects. Given that microRNAs have many more predicted targets, the number of proteins influenced by these 53 microRNAs is exponentially larger and has the potential to dramatically re-shape the living tissue environment.

Dlk1-Dio3 locus activation is characteristic of a subset of human lung adenocarcinoma samples

The small RNA population from 346 human lung adenocarcinoma samples has been subject to high throughput sequencing, available for analysis as part of the cancer genome atlas (TCGA). The percent of chr12qF1 microRNAs were increased 3.31-fold in these TCGA adenocarcinoma samples relative to 40 matched normal lung samples (p=0.014; Supplementary Figure 8a). However, from the tumor samples, a bimodal expression pattern was noted whereby several samples exhibited elevated expression patterns (Supplementary Figure 8a). Parsing of tumor samples into the 34 with high locus expression (~10%) with all other tumor samples revealed a consistent and specific activation of all Dlk1-Dio3 locus microRNAs (Supplementary Figure 8b). Thus, this locus is aberrantly activated in a subset of human lung adenocarcinomas.

Discussion

A central tenet of cancer genetics is that rapidly dividing tissues over time have the potential to accumulate enough mutations and chromosomal alterations in oncogenes and tumor suppressors such that a critical threshold is obtained and tumorigenesis ensues. However, in organs that do not undergo this rapid cell division and turnover – such as the lung – a reversal back to an embryonic state or a proliferation of stem cells (cancer stem cells) with their associated rapid growth and development is one of the mechanisms that is postulated to be involved. This is indicative of the scenario that we observe in this murine lung adenocarcinoma dataset in which a cluster of stem-cell associated microRNAs are up-regulated.

Several lines of evidence point to the role of the chr12qF1 locus in stem cell biology, lung development and oncogenesis. The sustained expression of this locus is essential in the development of induced pluripotent cells with a common loss of chr12qF1 microRNA expression resulting in the low proportion of cells that maintain an iPS state.[25] The proper expression of mmu-miR-127 is essential for lung development as its overexpression led to fewer terminal buds indicating impaired lung branching.[26] Meanwhile, removal of the maternally-derived (but not paternally-derived) Meg3 allele in mice led to thin-walled lungs with reduced radial alveolar counts and early postnatal lethality.[27] In addition, knockout of the Dicer1 gene, critical for RNA interference, had effects that were more specific to the lung,[28] while loss of one Dicer1 allele reduced survival in the same Kras lung adenocarcinoma model.[29] Several of the chr12qF1 microRNAs were among the most up-regulated by microarray analysis in a completely distinct mouse model of lung adenocarcinoma with sustained high levels of cyclin E and up-regulation of miR-376a and miR-136 was validated in human lung adenocarcinomas.[30] Differential expression of the chr12qF1 locus has also been identified in other cancer subtytpes[31] including up-regulation in mouse and human hepatocellular carcinoma samples,[32] gastrointestinal stromal tumors,[33] acute promyelocytic leukemia[34] and associated with epithelial to mesenchymal transition in endrometrial carcinoma.[35] Notably, miR-127-3p was significantly up in colorectal cancer associated with KRAS mutations.[36]

When an entire cluster of genes and miRNAs is up-regulated, it becomes difficult to identify and target a single gene or miRNA that may be responsible for the tumorigenic phenotype. It will be interesting to determine whether one or a few microRNAs are sufficient to recapitulate this oncogenic event, though a more likely scenario is that the cohort of small RNAs and non-coding RNAs act coordinately to regulate a multitude of genes. Several genes that are targets of chr12qF1 by luciferase analysis have implications in oncogenesis and lung development. The p53 interacting protein Tp53i11, repressed by miR-134 and miR-758 is a putative tumor suppressor in liver cancer.[37] The actin related protein Actr3 is a major constituent of the Arp2/3 protein complex down-regulated in gastric cancer,[38] Arfgef1 in breast cancer,[39] and Cul1 in various tissues.[40] Relevant to lung biology, the cystic fibrosis transmembrane receptor Cftr which is hypermethylated and down-regulated in lung adenocarcinoma[41] and mutated in non-small cell lung cancers,[42] is targeted by miR-381 and miR-494. Indeed several oncogenic pathways are implicated upon microRNA activation. How these microRNAs are specifically activated (or how a stem cell like population can continue to proliferate unabated) will be of interest in future studies. While methylation patterns do not appear to change in the tumor samples (contradictory to established methylation patterns at this locus), other epigenetic marks certainly may be involved, particularly histone marks. Alternatively since we examined methylation patterns of tumor samples in bulk, it remains possible that demethylation at a subset of cells is sufficient to activate the chr12qF1 locus. The results of this and other mouse models[43] appear to be limited to a subset of human lung tumors, and it will be interesting to determine if KRAS mutations in human samples can induce a similar activation of the chr12qF microRNAs.

By sequencing the transcriptome to a considerable depth we could search for additional mutations and expression changes that could aid in progression of tumorigenicity. An evaluation of both normal and tumorigenic lungs indicated no difference in the number of coding changes. However, if a variant caused nonsense-mediated decay or induces loss of expression, it would be more difficult to detect using this approach and would require whole genome or exome DNA sequencing. By examining allele-specific expression, we identified genes with modifications in expression of maternal versus paternal alleles depending on tumor status. This includes the paternal-specific enhancement of Cd22 mRNA expression specifically in tumor samples. CD22 is a B-cell lymphocyte cell surface marker that incurs mutations and splicing defects in human B-precursor leukemia[44] and is the target of Epratuzumab, a humanized monoclonal antibody therapeutic for B cell tumors.[45] Interestingly, a group recently independently identified CD22 cell surface expression in A549 cells and solid tumors and have shown that anti-CD22 antibodies can delay tumor progression.[46] Quantification of RNAseq read hits enabled a calculation of expression changes across the genome. The three most down-regulated genes (Prkg1, Gnao1 and Rgs9) are implicated in G-protein coupling and which have indirect connections with Kras and p53. Of note, a point mutation in Gnao1 has been recently identified in breast cancer,[11] and this mutation appears to function in a manner analogous to the Kras G12D mutation in that it maintains the gene in a constitutively active state.[47]

The presence of editing in these lung adenocarcinoma samples is quite striking, yet it is difficult to ensure that this effect is specific to tumor progression or infiltration of B cells in tumor samples. Of the 32 Apobec1 edited sites identified previously,[22] all but one are present in the 3′UTR of transcripts which can have multiple effects including influencing poly-A usage.[48] This can have regulatory consequences if regions such as microRNA binding sites are precluded from the edited RNA transcript. Our work suggests that much of the editing does not occur in tumor-originating cells; nonetheless, Apobec1 mRNA was increased overall in tumors and Apobec transgene expression in the liver has been shown to inadvertently drive HCC.[52] Separate from RNA editing, APOBEC mediated DNA mutagenesis from up-regulated APOBEC family proteins was reported to be a property of several human cancers[49, 50] after it was noted that cascades of localized C-to-T changes, termed kataegis, were found in breast cancer samples.[51]

In conclusion, through use of high throughput sequencing technology, we have uncovered several novel genetic abnormalities that exist in the coding and non-coding transcriptome of extracted solid tumors of the lung. By small RNA sequencing we were able to detect consistent up-regulation of a cluster of microRNAs typically associated with a stem cell like state. It is the activation of this locus and the multitude of mRNA targets of the ~53 microRNAs that we believe are crucial for oncogenic drive in this Kras mutant mouse model of lung adenocarcinoma.

Materials and Methods

Mouse breeding

The Stanford Institute of Medicine Animal Care and Use Committee approved all animal studies and procedures. 129S4 males heterozygous for a Kras allele were mated with MOLF/EiJ females. Littermate controls that were used include a Kras mouse with no Adeno-Cre administered, and a wildtype Kras mouse with Adeno-Cre (to control for adenovirus exposure effects). In addition, tumors and normal lungs were extracted from a distinct mouse lung adenocarcinoma model, namely the Kras LA mouse which does not require Cre delivery and is present on an inbred background.[8]

Kras mice were bred as previously described.[53, 54] To evaluate allele-specific expression, Kras mice on the 129S4/SvJae background were bred with MOLF mice. Offspring were genotyped for the presence of the Kras allele. To activate Kras, an adenovirus bearing Cre-recombinase (AdCre:CaPi co-precipitates, Baylor Vector Development Lab) was intra-nasally administered at 6 weeks of age. Tumors were dissected ~4 months later.

Tumor cell sorting and transplantation

The Rosa-LSL-tdTomato Cre reporter allele was bred into Kras p53flox mice. Tumors were initiated by intratracheal infection of mice with a lentiviral vector expressing Cre recombinase.[54] Single cell suspensions were generated from individual lung tumors harvested from mice 8–14 months after tumor initiation. Tumors were minced and then digested for 30 min at 37°C in 2ml of HBSS-free containing trypsin, collagenase IV and dispase. Subsequently, 4ml of Quench Solution (L15 media supplemented with FBS and DNase) was added and samples were then pressed through 40mm cell strainers (BD Biosciences). Finally, samples were centrifuged at 1,000 r.p.m. for 5 min and re-suspended in FACS media (PBS, 2% FBS, 2mM EDTA). Sorting for td-Tomato-positive cells was performed at the Stanford Shared FACS Facility.

High throughput large RNA sequencing

The Ribo-Zero ribosomal removal kit (Epicentre) was used to remove ribosomal sequences from total Trizol-extracted RNA. 200ng of purified RNA was subjected to strand-specific RNAseq using the ScriptSeq library preparation kit (Epicentre). 50bp paired end reads were generated on an Illumina HiSeq 2000 machine. Trimmed sequences were mapped to the mm9 genome using TopHat version 1.3.3 under default settings allowing multiple alignments with the following specific parameters: -r 100, --mate-std-dev 100 --segment-length 20–library-type fr-unstranded.[55] To accomodate SNPs between the Molf/EiJ and 129S4/SvJae strains, three mismatches were allowed in the mapping process for these samples, while two mismatches were allowed for samples with homozygous parents. FPKM values were calculated using Cufflinks v.1.3.0.[56, 57] Transcripts with a q-score of < 0.05 were considered significant and we further required a minimum FPKM of 5 in at least one of the two conditions.

SNP analysis

GATK (v.1.5) was used to calculate non-synonymous variants present in tumor and normal lung samples.[58, 59] Candidate variants in tumors were filtered based on a minimum read depth of 10, an absence in control lung samples, and that did not exhibit a bias in end-of-sequence. Raw reads were manually inspected for accuracy of GATK calls.

Allele specific expression

To evaluate allele specific expression, we extracted a list of all coding single nucleotide variants between MOLF and 129S4 mice as computed from the Mouse Genome Informatics dataset (Jackson labs). The Samtools BCFtools program[60] was used to call the number of variants in each sequenced sample. This data was filtered to have a minimum of 20 reads per sample. Variants in individual genes were sorted based on combined fold-changes between tumor and normal samples. Genes with the greatest fold-differences were manually analyzed for additional exonic SNPs.

RNA editing analysis

Allele calls at adenosine to inosine[21] and cytosine to uridine[22] editing sites were identified using BCFtools. The Serinc1 variant at chr10:57,235,791 (mm9) was validated by Sanger sequencing using the primers 5′-ACATTAGGCTCGGGTTAGGCACTA and 5′-AAGGCTGGAACATGAAGATGAACT for both genomic DNA and cDNA.

Small RNA sequencing

3ug of a mirVana (Invitrogen) extracted small RNA fraction was ligated to a 3′linker in ATP-free buffer. Samples were resolved on a 12% polyacrylamide gel and 17–28nt fragments were excised. A barcoded linker was ligated to the 5′ end of the extracted RNA using T4 RNA ligase and these RNAs were reverse transcribed using Superscript II (Invitrogen). 5ul of this product was subject to 21–24 PCR amplification cycles in a total of 50ul volume using Taq polymerase (NEB). The product was resolved on a 4% Nusieve GTG agarose gel (Lonza, Rockland, ME). 20ng of this product was subjected to 36 base pair high throughput sequencing on an Illumina GAII machine. Small RNA reads (17–28bp) were aligned to miRBase (release 15)[61] using the Bowtie program release 0.12.7[62] allowing for two mismatches.

Small RNA northern blot

10ug of Trizol (Invitrogen) extracted RNA was resolved on a 15% acrylamide gel, transferred to a Hybond N+ membrane (Amersham). Membranes were scanned using a phosphoimager.

Quantitative PCR

2ug of total RNA was reverse-transcribed using Superscript II (Invitrogen). Gene-specific probes used included Rtl1 (Mm02392620_s1), Dlk1 (Mm00494477_m1), Meg3 (Mm00522599_m1), Tnrc6b (Mm00523487_m1), Cftr (Mm00445197_m1) and Antxr1 (Mm00712952_m1) (Life Technologies). MicroRNA Taqman analysis was performed on 250ng of Trizol-extracted small RNAs using probes that detect mature miR-21, miR-127-5p, miR-376a and U6 (Life Technologies). Quantitative RT-PCR was performed on a CFX384 Real-Time system (BioRad).

Luciferase analysis

Dual-luciferase assays (Promega) were performed 24hr after transfection according to manufacturer’s protocol and detected by a Modulas Microplate Luminometer (Turner Biosystems). For transfection, 250ng of psi-check reporter plasmids were co-transfected with 250ng of miRNA overexpression plasmids (Sh-constructs) in E10.5 mouse embryonic fibroblasts using TransIT-LT1 (Mirus). Cell seeding was performed at a concentration of 2.5×104 cells per well in a 24-well plate.

For cloning of the psi-check constucts, the entire 3′UTR of each gene was PCR amplified (see Supplementary Table 5 for primer sequences) from mouse genomic DNA and cloned in psi-check-2 vector (Promega) between the XhoI and SpeI sites using the In-fusion HD cloning kit. The quickchange II site directed mutagenesis kit (Agilent) was used to introduce three mutations in the Hyal1 3′UTR using the primer 5′-GGACTTCCTCAAATACTGACTCATGCCCATAAGTC and the reverse complement thereof (mismatches are listed in bold). For generation of the microRNA over-expression constructs, shRNA sequences (Supplementary Table 6) were chemically synthesized; both strands were annealed and inserted between BglII and KpnI sites downstream of the U6 Pol III promoter.

Bisulfite sequencing

Genomic DNA was extracted from tumor and normal lung samples by the DNeasy tissue kit (Qiagen) and converted by bisulfite treatment using an EZ DNA methylation kit (Zymo research). PCR primers were as follows for Meg3: 5′-GTTATAGTAATTTGTTATAGAATTTGGGG (forward) and 5′-AAACTTTCAACCACCAAAACC (reverse), and for an intergenic differentially methylated region: 5′-GGTTTGGTATATATGGATGTATTGTAATATAGG (forward) and 5′-ATAAAACACCAAATCTATACCAAAATATACC (reverse).[63] Products were cloned using the Topo TA cloning system (Invitrogen) and at least eight clones were sequenced per sample per locus.

miRNA target prediction

For Miranda, 395 targets were identified for chr12qF1 microRNAs versus 376 in controls using a cutoff mirSVR score of -1.3[64]. The Pictar program predicted 1231 chr12qF1 and 1171 control target mRNAs. Finally, 1002 chr12qF1 and 1125 control mRNAs were identified as TargetScan [61, 65] targets (within the 95th percentile of hits).

Pathway analysis

Transcripts significantly differentially expressed between tumor and normal samples (false discovery rate < 0.001) were analyzed by the ingenuity pathway analysis program using default settings.

Human lung adenocarcinoma samples

Small RNA deep sequencing reads corresponding to human microRNAs were queried from samples collected as part of the cancer genome atlas (TCGA). Data was available from a total of 346 samples from patients with stage I–IV lung adenocarcinomas collected from centers throughout the United States. This was compared with data from 40 matched normal samples.

Data access

Sequences have been deposited in the NCBI Gene Expression Omnibus (accession number GSE43028).