Deregulated lncRNA expression profile in the mouse lung adenocarcinomas with KRAS‐G12D mutation and P53 knockout

Abstract Recent studies have demonstrated that aberrant long non‐coding RNAs (lncRNAs) expression are suggested to be closely associated with multiple human diseases, lung cancer included. However, the roles of lncRNAs in lung cancer are not well understood. In this study, we used microarrays to investigate the aberrantly expressed lncRNAs in the mouse lung adenocarcinoma with P53 knockout and the KrasG12D mutation. Results revealed that 6424 lncRNAs were differentially expressed (≥ 2‐fold change, P < .05). Two hundred and ten lncRNAs showed more than 8‐fold change and conserved across human and were further analysed in the primary mouse lung adenocarcinoma KP cells, which were isolated from the p53 knockout and the KrasG12D mutation mice. Among all the 210 lncRNAs, 11 lncRNAs' expression was regulated by P53, 33 lncRNAs by KRAS and 13 lncRNAs by hypoxia in the primary KP cells, respectively. NONMMUT015812, which was remarkably up‐regulated in the mouse lung adenocarcinoma and negatively regulated by the P53 re‐expression, was detected to analyse its cellular function. Results showed that knockdown of NONMMUT015812 by shRNAs decreased proliferation and migration abilities of KP cells. Among those aberrantly expressed lncRNAs in the mouse lung adenocarcinoma, NONMMUT015812 was a potential oncogene.


| INTRODUC TI ON
As the most prevalent primary malignant tumour in the world, lung cancer has been the leading cause of cancer-related death. 1 Nonsmall cell lung cancer (NSCLC) is the most common type of lung cancer and includes adenocarcinoma, squamous cell carcinoma and large cell carcinoma. 2 Adenocarcinoma is the most common type of NSCLC, now accounting for about 40% of all lung cancer cases. 3 Genetic alterations in non-small cell lung cancer (NSCLC) have been identified in recent years. Kristen rat sarcoma viral oncogene (KRAS), epidermal growth factor receptor (EGFR) and anaplastic lymphoma kinase (ALK) are the most commonly altered oncogenes by acting as tumour genomic drivers. 4 More recently, the identification of activating epidermal growth factor receptor (EGFR) mutations and anaplastic lymphoma kinase (ALK) rearrangements as predictive biomarkers for treatment of NSCLC led to further personalization of therapy with EGFR tyrosine kinase inhibitors (EGFR-TKIs) and ALK inhibitors. KRAS is mutated to an activated form in ∼30% of NSCLCs; however, effective clinical therapies targeting KRAS have yet to be developed. The study of the underlying biology of KRAS in patients with NSCLC could help to determine potential candidates to evaluate novel targeted agents and combinations that may allow a tailored treatment for these patients. 5 In addition, the P53 tumour suppressor gene is mutated or deleted in ∼50% human cancers. 6 Wild-type (WT) P53 helps maintain genome integrity and cellular homeostasis by regulating the expression of a plethora of genes involved in the regulation of cell cycle, apoptosis, stem cell differentiation, senescence, DNA repair and metabolism. 7 Long non-coding RNAs (lncRNAs), a kind of transcript RNA molecular longer than 200 nucleotides, are not translated into protein in the nucleus or cytoplasm. Studies have shown that lncRNAs play significant roles in many life activities such as epigenetic regulation, cell cycle regulation and cell differentiation regulation, and abnormally express in a variety of cancers. The aberrant expression of lncRNAs can be used as tumour promoters or inhibitors. Studies have indicated that lncRNAs play critical roles in lung cancer development and progression by altering multiple signalling pathways, and some cancer-associated lncRNAs have been identified, such as H19, 8 MALAT(lung adenocarcinoma associated transcript 1), 9,10 MIR31HG, 11 HOTAIR, 12 DGCR5, 13 AFAP1-AS1, 14 SNHG1 and RMRP. 15 However, the research on the identification of additional lung cancer-associated lncRNAs remains to be investigated.
The widely used lung adenocarcinoma mouse model with p53 deletion and Kras G12D mutation was chosen in this study. The expression patterns of lncRNAs in the mouse lung adenocarcinoma tissue were detected by using microarrays. Based on the microarray results, 210 conserved lncRNAs were selected to detect the relationship between their expression and P53 or KRAS. Besides, hypoxia is an important tumour microenvironmental factor that can induce and even aggravate metastasis, 16 so we also detected the effects of hypoxia on the 210 lncRNAs. Furthermore, the cellular function of NONMMUT015812 was further investigated in the primary mouse lung adenocarcinoma cells.

| Mouse tissue samples
The Institutional Animal Care and Use Committee of Fudan University, China, approved all protocols. p53(flox/flox) mice (The Jackson Laboratory) were crossed with the LSL-Kras G12D mice (The Jackson Laboratory) to generate p53(flox/flox); LSL-Kras-G12D mice and the genotyping primers were listed in the Table S1.
As described by DuPage, 17 lung adenocarcinomas were initiated by intranasal instillation with 1 × 10 6 lentiviral particles expressing Cre-recombinase per 8-week-old mouse. The lung adenocarcinoma samples and corresponding normal lung samples were prospectively collected 6 weeks after tumour initiation under a dissection microscope (Table S2). The tumour size ranges from 2 to 4 mm and, consistent with DuPage's description, these tumours were classified as Grade 2. No obvious metastatic lesions were found in livers and brains.

| Microarray and computation
For microarray analysis, lncRNA + mRNA mouse Gene Expression Microarray v1.0, 4x180K chip (CapitalBio Corp), was employed and conducted by CapitalBio Corp according to the manufacturer's protocols. Briefly, 1 μg of total RNA extracted from the samples was transcribed into double-stranded cDNA using CbcScript reverse transcriptase and T7-random primer according to the manufacturer's protocol (CapitalBio). Double-stranded cDNA products were purified and in vitro transcribed to cRNA with T7 RNA polymerase.
After reverse transcription, the cRNA was transcribed into cDNA and labelled with Cy3 by using the Klenow enzyme-labelling strategy, and microarray slide hybridized with the Cy3 labelled probes.
Feature Extraction Software v10.7 was used for data extraction from raw microarray images files. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v13.0 software package (Agilent Technologies).

| RNA extraction and quantitative real-time PCR
Total RNA was extracted from cultured cells or tissues using TRIzol reagent (Thermo Fisher) according to the manufacturer's protocol.
RNA degradation and contamination were assessed on 1% agarose gels, and RNA concentration was measured using a NanoDrop 1000 (Thermo Scientific). cDNA was synthesized using PrimeScript RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa), and the integrity of synthesized cDNA was confirmed using glyceraldehyde 3-phosphate dehydrogenase (Gapdh) or 18s as the endogenous control. Real-time PCR was carried out using SYBR Premix Ex TaqTM II (Perfect Real Time) (TaKaRa) and measured using ABI 7900 or ABI 7500 instrument. The primers are listed in Table S3. PCR was per-  The Kras CRISPR guide RNA sequences were designed by using the web tool developed by Feng Zhang group at MIT and subcloned between the two BsmBI sites of the Lenti-Guide-coGFP-P2A-Puro guide RNA Expression Lenti-vector. Plasmid lentiCas9-Blast from Dr Feng Zhang's laboratory was used to package Cas-expression lentivirus, and PCDH-Cre was used to package Cre-expression lentivirus.

| Plasmid expression constructs
All the used DNA oligos were listed in the Table S1. The infected stable KP cells were selected with 4 µg/mL puromycin dihydrochloride (Amresco) or 10 µg/mL blasticidin (Invitrogen) after 72 hours by using lentiviral infection.

| RNA sequencing
The NONMMUT015812-knockdown KP (shRNA-2) cells and negative control (sh-Scr) cells were obtained on day 6 of a viral infection. In-depth RNA sequencing was performed by the Biomarker Technologies in China. The raw sequencing image data were examined via Illumina analysis pipeline and aligned with the unmasked mouse reference genome. Differential expression analysis of two groups was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < .01 and absolute value of log2(Fold change)>1 found by DESeq were assigned as differentially expressed.

| Construction of lncRNA-lncRNA co-expression network and statistical analyses
Based on the computational methods used by Pan JQ, 18 the expression values of differential expressed lncRNAs (DE lncRNAs) were obtained from Table S4. The correlation coefficient (CC) and significant P-value were calculated between the expression values of each lncRNA-lncRNA pair across all samples. The lncRNA-lncRNA pairs with a P-value < .001 and absolute value of CC > 0.9 were used to construct a lncRNA-lncRNA co-expression network (LLCN). The LLCN was generated by using Cytoscape software.
Statistical differences (P-values) among groups were obtained using two-sided Student's t test. Results were considered to be significant for P values under .05. All statistical tests were performed using Prism software (GraphPad, version 5.0). The Kras mutation gene (c.35G > A) was confirmed by DNA sequencing ( Figure 2B). Western blot showed that P53 was completely deleted in the KP cells ( Figure 2C). The primary KP cells derived from p53(flox/flox);LSL-Kras-G12D C57 mice, so they could theoretically form tumour in wild-type C57. To assess the capability of tumour formation, KP cells were subcutaneously injected into armpits of wild-type C57 mice, and obvious tumour formed within one month and adenocarcinoma was confirmed by the H&E staining analysis ( Figure 2D and 2E). Among the selected 210 lncRNAs, qPCR results determined 11 lncRNAs, whose expression was altered by the re-expression of P53 in KP cells and exhibited opposite expression alteration direction to that in the microarray of the mouse lung adenocarcinoma tissues. The heights of the columns in the chart represent the log-transformed median fold changes (tumour sample/normal tissues in microarray data or KP-p53 cells/KP-NC cells in qPCR results). Mean ± SD are shown, n = 3 or 2

| Analysis of lncRNAs associated with p53 Kras or hypoxia
p53 is the most common tumour suppressor gene mutated in clinical lung cancer, and it was induced to complete knockout in our lung adenocarcinoma mouse model and the KP cells. Thus, we explored whether the deregulated expression of the selected conserved 210 lncRNAs was associated with the p53 deletion. P53 was re-expressed in the KP cells by the lentiviral particles carrying the p53 gene. As shown in Figure 3A, P53 was efficiently re-expressed in the stably infected KP cells. The quantitative PCR assay showed that expression of 11 lncRNAs from the selected 210 lncRNAs was significantly altered by P53 re-expression (fold change ≥ 1.5), and the fold change was contrary to the above microarray data of lung adenocarcinoma tissues in which p53 was deleted ( Figure 3B).
Kras was activated in the lung adenocarcinoma mouse model and the KP cells. Similarly, to identify lncRNAs whose expression was relevant to the Kras status, we used the Crispr/Cas9 targeting the coding sequence region of Kras gene to knockout its expression in KP cells.
The DNA sequencing of genomic PCR products containing the guider RNA binding site showed that the Kras gene was efficiently mutated in KP cells by the Crispr/Cas9, and the cell growth was obviously inhibited ( Figure 4A). The quantitative PCR assay identified 33 lncRNAs whose expression was affected by the deletion of Kras (fold change ≥ 1.5), and the fold change was also contrary to the original microarray data of lung adenocarcinoma tissues in which Kras was activated ( Figure 4B).

| Knockdown of NONMMUT015812 suppressed KP cell proliferation and migration
Among the identified lncRNAs, NONMMUT015812 was very interesting. (a) Its expression increased approximately 150-187-fold in the mouse lung adenocarcinoma tissues. (b) Its expression was negatively associated with P53. (c) The gene for NONMMUT015812 is located on the chromosome between the genes for retrotransposon-like protein 1 and thyroxine 5-deiodinase, in which several differently expressed lncRNAs were included and conserved in human and mouse (this region will be described in detail in the discussion). There were no reports about NONMMUT015812, so we further studied its cellular function in the KP cells. Lentivirusmediated shRNAs were designed to knockdown the expression of NONMMUT015812, and the results of quantitative PCR showed that they could reduce NONMMUT015812 in KP cells by more than 80% ( Figure 7A). Clone forming assays showed that knockdown of NONMMUT015812 significantly inhibited KP cell growth ( Figure 7B). Moreover, in transwell assays, KP cells with stable NONMMUT015812 knockdown showed significantly decreased migration compared with control cells (Figure 7C). Flow cytometry F I G U R E 5 Identification of hypoxia-related lncRNAs in KP cells. A, KP cells were incubated under hypoxia with 1% O 2 and normoxia (20% O 2 ). Immunoblot (left) was used to detect HIF1α protein levels and was quantified by measuring the grey value of the bands. Data are shown as means ± SD. All P values are from two-sided t tests. B, Comparison between microarray data and qPCR results of lncRNAs associated with hypoxia. The qPCR result validated the 13 out of 210 selected lnRNAs were regulated by hypoxia, whose expression exhibited same expression alteration direction to that in the microarray of the mouse lung adenocarcinoma tissues. The heights of the columns in the chart represent the log-transformed median fold changes (tumour sample/normal tissues in microarray data of the mouse lung adenocarcinoma tissues or KP cells cultured in hypoxia/KP cells cultured in normoxia in qPCR results) in expression across the three lung adenocarcinoma tissues and KP cells cultured in hypoxia. Mean ± SD are shown, n = 3 or 2 analysis indicated that NONMMUT015812 knockdown resulted in a substantial retardation of KP cells in G0/G1 phase, with decrease in the number of cells in G2/M phase ( Figure 7D). However, knockdown of NONMMUT015812 had no significant effects on the apoptosis ( Figure S1).
RNA sequencing was used to monitor the change in gene expression in the NONMMUT015812-knockdown cells. Compared with the negative control, 1281 differentially expressed (DE) genes were identified. Among those genes, 610 genes were up-regulated, and 671 genes were down-regulated ( Figure 7E, Table S8 and S9). KEGG pathway analysis indicated that down-regulation of NONMMUT015812 affected several cell growth-associated pathways, including cell cycle, DNA replication and PI3K/Akt signalling pathway ( Figure 7F, Table S10). Interestingly, P53 signalling and several DNA repair pathways were also included in the KEGG pathway enrichment, which means that NONMMUT015812 might be an important downstream mediator of P53 pathway. Consistent with the results of flow cytometry analysis, no apoptosis signalling pathway was identified. Gene Ontology analysis also revealed that many differentially expressed genes were mainly enriched in cell division, mitotic cell cycle and DNA replication processes ( Figure 7G, Table S11).
F I G U R E 6 Confirmation of the p53-, Kras-or hypoxia-associated lncRNAs. A, Overlapped differentially expressed lncRNAs. B, Among the total 52 lncRNAs identified to be associated with p53, Kras or hypoxia, 19 lncRNAs, with expression change fold greater than 3.5 in the microarray results, were selected to further confirmed by using quantitative RT-PCR in another group of mouse lung adenocarcinoma tissues (6 lung adenocarcinoma samples and 6 normal lung tissues). The heights of the columns in the chart represent the log-transformed median fold changes (tumour sample/normal tissues). The validation results of the 19 lncRNAs indicated that the microarray data correlated well with the qPCR results. Mean ± SD are shown, n = 3 or 6. C, Long noncoding RNA (lncRNA)-lncRNA co-expression network (LLCN) related to the 19 selected lncRNAs was constructed. Among the 19 lncRNAs, 6 lncRNAs had no significant co-expression with other lncRNAs. The circles represent nodes of differential expressed lncRNAs (DE lncRNAs). The nodes in green colour represent the 13 remained lncRNAs. The LLCN contained a total of 4086 lncRNAs with 19 595 edges between them

| D ISCUSS I ON
Lung cancer is the leading cause of cancer death, and the majority of cases are diagnosed as NSCLC, in which lung adenocarcinoma is the most common type. 23 Recently, significant progress has been made in individualizing therapy based on molecular aberrations and pathologic subtype. For example, the use of EGFR-tyrosine kinase inhibitors or ALK inhibitors has greatly improved progression-free survival in NSCL patients with EGFR mutations or ALK translocations. KRAS and P53 are among the most frequently mutated genes in non-small cell lung cancer, but no targeted therapy has been approved for KRAS mutant or P53 mutant NSCLC at present. Therefore, a better understanding of the pathogenesis of lung adenocarcinoma with KRAS or P53 mutation could provide more effective diagnostic markers or treatment targets for lung adenocarcinoma.
Recent investigations have provided good evidence that lncRNAs play important roles in many biological processes, from chromatin modification, transcription, splicing and translation to cellular differentiation, cell cycle regulation and stem cells reprogramming. 24 Furthermore, some lncRNAs have been demonstrated to be aberrantly expressed in lung cancer tissues and be associated with cancer progression. For example, HOTAIR is up-regulated in lung cancer specimens and displayed potential tumour-promoting roles 25-27 ; MALAT1 is a critical regulator of the metastasis phenotype of lung cancer cells, as a competing endogenous RNA (ceRNA) to regulate autophagy-related 7 (ATG7) gene expression by sponging miR142-3p. 28 To investigate lncRNAs deregulated in p53-and Kras-mutated lung cancer, the microarrays were used to analyse lncRNA expression profiles in the tissues of mice lung adenocarcinoma tissues with Besides, some lncRNAs showed to be regulated by two factors.
We noticed an intriguing phenomenon that 16 of the identified conserved lncRNAs located on the chromosome region between the genes for retrotransposon-like protein 1 and thyroxine 5-deiodinase, they were all up-regulated in the mouse lung adenocarcinoma tissues, and expression of 9 lncRNAs was identified to be associated with p53, Kras mutation or hypoxia. The lncRNA maternally expressed 3 (MEG3) reported to be correlated with tumorigenesis, 29,30 was also included in this region and was identified to be up-regulated in the mouse lung adenocarcinoma tissues (corresponding to the lncRNA NONMMUT015697) and regulated by kras in our results. The above results suggested that this region contain a cluster of tumorigenesis-promoting lncRNAs.
The lncRNA NONMMUT015812, which our results showed to be highly up-regulated in the mouse lung adenocarcinoma by more than 150-fold and suppressed by P53, was another lncRNA in this region, and there was no research on its function. Our functional analysis showed that knockdown of NONMMUT015812 suppressed KP cell proliferation and migration, which means it was a candidate oncogene. The RNA-seq further indicated that NONMMUT015812 mainly regulated cell cycle to control cell growth( Figure 7).
In conclusion, the deregulated lncRNA expression profile was determined in the mouse lung adenocarcinomas with Kras-G12D mutation and p53 knockout. The expression of some highly dysregulated and conserved lncRNAs was revealed to be associated with p53, Kras or hypoxia, among which the lncRNA NONMMUT015812 was a potential tumorigenesis-promoting gene. However, the function and molecular mechanism of the identified lncRNAs in tumorigenesis and their expression in human lung cancer tissues need further investigating.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest.
F I G U R E 7 NONMMUT015812 is a potential oncogene. A, NONMMUT015812 was knocked down in the KP cells. Relative expression levels of NONMMUT015812 were determined by qRT-PCR in the KP cells infected with lentivirus of shRNA-1 and shRNA-2 to target NONMMUT015812, with the non-specific shRNA as negative control. Results are shown as means ± SD. B, Knockdown of NONMMUT015812 inhibited colony formation in KP cells. C, Transwell migration assay of the KP cells stably expressing shRNAs. Compared with the non-specific shRNA as negative control, knockdown of NONMMUT015812 significantly inhibited the cell migration. D, Flow cytometry assays showed that NONMMUT015812 knockdown inhibited the cell cycle progression in KP cells. E, A volcano plot of differentially expressed genes in the NONMMUT015812-knockdown KP cells. RNA sequencing revealed that a total of 1281 genes were differentially expressed between the NONMMUT015812-knockdown KP cells and the negative control (up-regulated: 610, down-regulated: 671). F, KEGG pathway analysis of the differentially expressed genes in the NONMMUT015812-knockdown KP cells. G, Gene Ontology analysis of the differentially expressed genes in the NONMMUT015812-knockdown KP cells, which revealed that many differentially expressed genes were mainly enriched in cell division, mitotic cell cycle and DNA replication processes (P < .001)

AUTH O R CO NTR I B UTI O N S
Jin Zhang and Duan Ma designed the study. Meiqin Zhang performed the experiments and was primarily responsible for writing the manuscript. All the authors contributed to data analysis and manuscript editing and approval.