Experimental Verification of a Predicted Intronic MicroRNA in Human NGFR Gene with a Potential Pro-Apoptotic Function

Neurotrophins (NTs) are a family of secreted growth factor proteins primarily involved in the regulation of survival and appropriate development of neural cells, functioning by binding to their specific (TrkA, TtkB, and TrkC) and/or common NGFR receptor. NGFR is the common receptor of NTs, binding with low-affinity to all members of the family. Among different functions assigned to NGFR, it is also involved in apoptosis induction and tumorigenesis processes. Interestingly, some of the functions of NGFR appear to be ligand-independent, suggesting a probable involvement of non-coding RNA residing within the sequence of the gene. Here, we are reporting the existence of a conserved putative microRNA, named Hsa-mir-6165 [EBI accession#: FR873488]. Transfection of a DNA segment corresponding to the pre-mir-6165 sequence in Hela cell line caused the generation of mature exogenous mir-6165 (a ∼200,000 fold overexpression). Furthermore, using specific primers, we succeeded to detect the endogenous expression of mir-6165 in several glioma cell lines and glioma primary tumors known to express NGFR. Similar to the pro-apoptotic role of NGFR in some cell types, overexpression of pre-mir-6165 in U87 cell line resulted in an elevated rate of apoptosis. Moreover, coordinated with the increased level of mir-6165 in the transfected U87 cell line, two of its predicted target genes (Pkd1 and DAGLA) were significantly down-regulated. The latter findings suggest that some of the previously attributed functions of NGFR could be explained indirectly by co-transcription of mir-6165 in the cells.

There are some transcription factors which regulate NGFR gene expression in conditions like hypo-osmolar stress and injuries [16]. There are also an increasing list of ligands, co-receptors and adaptor proteins which interact with NGFR, involving this gene in different signaling pathways ending to cell death and/or survival [17]. Recent publications suggest a ligand independent activation and also a non-NGF ligand activation of NGFR signaling [18][19]. It remains to be found if non-coding RNAs such as miRNAs are involved in this regulation as well.
miRNAs are endogenous small non-coding RNAs about 21-23nt, each capable of interfering with dozens of target mRNAs through complete or partial complementarities. The miRNA genes are transcribed by RNA polymerase II or III and the primary transcript termed pri-miRNA, is quickly trimmed into ,70nt long pre-miRNA precursor. The pre-miRNAs are then transported to the cytoplasm and further processed to functional miRNAs, which exert their regulatory functions by complementary binding to their target genes [20]. To date, more than 1000 human miRNAs have been published in mirbase database [21], among which 40% are located within introns of protein-coding 'host' genes [22][23][24][25][26][27] and appear to be conserved across the species.
Identification of novel miRNA by forward genetics has been tedious as a result of the small sizes of miRNAs as well as their tolerance to mutations that do not affect their seed sequences [28]. On the other hand, the computational prediction of non-coding RNAs (ncRNAs) has proven to be fast, cheap and effective [29][30]. After prediction of a miRNA gene, experimental verification is necessary to demonstrate its exact mature sequence and function. Briefly, some major characteristics such as the hairpinshaped stem loop structure, high evolutionary conservation and high minimal folding free energy are important features used in the computational identification of novel miRNAs [31]. Here bioinformatics tools were used to search for hairpin structures within the NGFR intronic and exonic regions. One of the NGFR resided hairpin structures showed all the bioinformatics characters of producing a real miRNA. Later, accumulative experimental evidences for the first time showed the presence and the functionality of this novel NGFR intronic miRNA in human cells.

Tissue Samples
Fresh surgical tissue biopsies of meningioma, glial, astrocytoma and oligodendroglioma were kindly provided by Imam Hospital [49].The samples were stored and processed as previously described in the same reference.

RNA Extraction
Total RNA was extracted from cell lines using Trizol reagent according to the manufacturer's protocol (Invitrogen). Residual DNA was removed using RNAase-free DNAase I (Takara) at 37uC for 30 min followed by heat inactivation at 65uC for 10 min by addition of EDTA.
Real-time PCR detection of precursor and mir-6165 mature form.
According to the predicted precursor and mature miRNA sequences, primers were designed for quantitative PCR (Table 1), using NCBI primer-blast, MWG Operon online PCR primer design tool (www.eurofinsdna.com) and primer bank (http://pga. mgh.harvard.edu/primerbank/). Briefly, 1 ug of total RNA Figure 2. pre-mir-6165 overexpression in the Hela cells and detection of Hsa-mir-6165 mature form. A) Schematic presentation of premir-6165 cloning, overexpression and its RNA adenylation followed by cDNA synthesis, using a universal anchored-oligo-dT primer. For the amplification of precursor, first strand cDNA was PCR amplified using precursor specific F-primer and reverse anchor primer on the oligo-dT tail. For the amplification of mature miRNA, predicted mir-6165 sequence was used as the forward primer. B) Hsa-mir-6165 increased production (200,000x) following transfection of Hela cells with its precursor. In the untransfected cells (U) or scrambled control (M), the level of this miRNA was lower compared to the transfected cells (T). C) Four sequencing result of TA vector clones containing mir-6165 real time PCR products, are compared to the precursor sequence. The sequences between the laboratory added polyA and the upstream vector sequence are considered as mature miRNA. Sequencing of clones #2, #3 and #4 shows that prediction of Hsa-mir-6165 sequence has been correct. Clone#1 also shows the similar sequence plus AGG extra nucleotides which is considered an iso-mir for the miRNA. doi:10.1371/journal.pone.0035561.g002 sample was used in 20 ul Poly A tailing reaction containing; 2.5 U PolyA polymerase (BioLab), 2 ul of 10 mMol ATP and incubated at 37uC for 10 min. Then, the total poly adenylated RNA was used in 10 mL first strand cDNA synthesis reaction using PrimeScript II reverse transcriptase (Takara) and a cocktail of specific oligo-dT primers containing a universal anchor (Table-1). cDNA synthesis reaction was performed at 42uC for 30 min and terminated at 80uC for 5 s. Real-time quantitative PCR was performed using standard protocols on an ABI PRISM 7500 instrument (Applied Biosystems). Briefly, the run method profile consisted of: stage 1, 95uC for 5 s,stage 2, 60uC for 20 s; stage 3, 72uC for 34 s. Stage 2 was repeated for 45 cycles. Continuous melt curve stages included a first step of 95uC/15 s, step 2 at 60.0uC for 1 min, 95uC/30 s and a last step of 60uC for 15 s. Total PCR products were cloned in TA vector (Fermentas) and were sequenced (Genfanavaran Co.).

Overexpression of mir-6165 Precursor in Hela Cell Line
Human genomic DNA was extracted from white blood cells using standard protocol [50]. 84 bp fragment representing the predicted pre-miRNA was PCR amplified, using precursor forward and reverse primers (Table 1), using PFU polymerase (Takara) and cloned in modified pEGFP-C1 expression vector with two CMV promoters. The 84 bp fragment was cloned in the vector using Not1 and EcoRV restriction enzymes. Recombinant vectors were propagated by miniprep (Qiagene Co.) and 2 ug of this DNA was used for lipofectamin (Invitrogen) transfection of Hela and U87 cell lines in 24 well plates containing about 26105 cells per well. GFP expression was visualized by a florescence microscope (Nikon eclipse Te2000-s).

Cell Cycle Analysis
After 34 hours, the U87 and Hela cells over expressing pre-mir-6165 were harvested and stained with propidium iodide (PI) and Annexin V (Roche) according to the manufacturer's protocol. All samples were analyzed with a FACS Calibur flow cytometer with Cell Quest software (BD Biosciences). All assays were carried out in duplicates.

Statistical Analysis
Real time experiments were run in duplicates. Real time data were analyzed using DDCT method by DataAssist software V3.0 and normalized by endogenous control U48 small nucleolar RNA gene (SNORD48) and B2m or globalization method [51]. Other statistical analysis was performed with GraphPad Prism 5.04 (GraphPad, San Diego, CA). For apoptosis studies, data showing percent of early apoptotic cell population within negative group and test group, compared with each other by Repeated Measures ANOVA test, and followed by Bonferroni test using GraphPad. Data were considered statistically significant, when P-values were ,0.05.

Results and Discussion
Prediction of a Novel Intronic miRNA Within the Human NGFR Gene Due to their implications in several diseases including cancer, miRNAs are under intensive research aiming at novel pharmacological interventions. Computational tools have been used for efficient prediction of novel miRNAs and their target genes [30]. After several failed PCR attempts to amplify the area spanning the 4 th intron of rat NGFR gene, chr10: 84267661-84265746[-] UCSC nov.2004 (Baylor 34 rn4), we used mFOLD program (http://mfold.rna.albany.edu) to search for possible hairpin structures in this region. This program introduced multiple hairpin structures in the 4 th intron of NGFR genes both in human and rat genomes. Using miRNA prediction HMM based tool (SCC profiler), several potential miRNA precursors were identified, however, only one of them, named pre-mir-6165, showed the criteria of producing a real human intronic miRNA locating; hg19, chr17: 47588166 -47588269 ( Figure 1A). Microprocessor SVM program also predicted a Dorsha processing site for pre-mir-6165 ( Figure 1B). To date for mir-6165, no identical miRNA has been reported in the mirbase database, except a weak resemblance to mir-328a in rat and mouse. Furthermore, Mireval online tool identified this novel precursor with strong conservation and with no homology to other miRNAs. Using blast search, it was demonstrated that both sequence and structure of pre-mir-6165 was conserved in mammals, while mir-6165 was mostly conserved in primates ( Figure 1C). Overall, accumulated bioinformatics evidences suggested the existence of a novel microRNA. Firstly, CIDMIR, Pmirp, MatureBayes, miRabela -MirZ and Microprocessor SVM softwares all recognized pre-mir-6165 with significant scores. These tools consider conservation, expression level and other characters in order to predict the possibility of a miRNA production. Secondly, precursor sequence as a query produced 9 hits with 100% identity at the level of 20-27 bp in the Blat search against human genome. Four of these yet uncharacterized hits show common characters of mir-precursors and are very similar to the already reported miRNA precursors in the mirbase database. Others have used similar method for the discovery of new miRNAs [32]. Thirdly, conservation of seed [52] as well as the rest of mir-6165 sequence in human genome is a strong bioinformatics supporting evidence for the presence of this miRNA ( Figure 1C). This miRNA is not clustered and is weakly conserved between taxa, the same property is already reported for certain miRNAs which are not conserved [53], not clustered [54], have temporal and cell-specific expression patterns [55], and have no homology to other miRNAs [35].

Overexpressed pre-mir-6165 is Efficiently Processed in the Hela Cells
According to the Oncomine [56], the database for the expression profile of the genes, NGFR gene expression (and as a result pre-mir-6165) is very low in Hela cell line. In an attempt to overexpress pre-mir-6165 in Hela cells and check for the production of mature mir-6165, a corresponding 84 bp DNA fragment was amplified by PCR from human genomic DNA and ligated to the pEGFP-C1 modified expression vector under CMV promoter (Figure 2A). The same vector without insert and/or containing scramble DNA, were used as a negative controls in further experiments. Transfection rate of these cells was monitored via fluorescent emission of GFP, and the best transfected culture was used for RNA extraction and gene expression studies. In Hela cells which were transfected with the pre-mir-6165 construct the level of mature mir-6165 elevated by a factor of 200,000 fold. The result indicates that mir-6165 precursor cloned under CMV promoter is efficiently expressed and processed to a mature form ( Figure 2B). Mature mir-6165 real-time PCR products showed an expected size on the acrylamid gel and were subsequently cloned in the TA vector for sequencing. The sequencing results of three independent clones showed the exact sequence of expected mir-6165 before polyA tail, indicating correct prediction of 39-end of the miRNA ( Figure 2C). On the other hand, sequencing result from one of the TA clones showed 2 bp extra nucleotides after the predicted 39 end for mir-6165. This might be the result of subtle variation in the RNA ends corresponding to Drosha and Dicer cleavage sites.
The minimum size of these miRNAs was submitted to EBI data base under the extension number of; [EBI accession#: FR873488, FR873489]. Overall, these experiments demonstrate that Hela cells are able to process the predicted pre-mir-6165 to its predicted mature miRNA form.

Experimental Detection of Endogenous mir-6165 and its Precursor in Glioma Cell Lines and Brain Tumors
Promoter 2.0 Prediction Server did not predict an independent promoter for mir-6165, therefore, its transcription is supposed to be through NGFR host gene promoter. NGFR is expressed shortly   in few normal cell lineages during the development, but it is expressed at a detectable level in the human Glioblastoma U87-MG cell line [15,56]. For this reason, U87-MG and other glioblastoma cell lines were chosen for the experimental detection of mir-6165 and its precursor. The fact that one or both of mir and miRNA-star sequence (miRNA*), the complementary strands of functional mature miRNA, might exert function, two primers were designed; one exactly identical to the predicted mir-6165 (called Pm1 primer, Table-1) and the second primer was identical to its complementary mir* (called Pm1* primer). After Poly-A addition to the total RNA extracted from several glioblastoma cell lines, cDNA was synthesized using anchored oligodT primer (Table-1). Later along with NGFR gene, the mature form of mir-6165 was amplified from RNA samples of Daoy, 1321N1, U-87 and NT2 cell lines, using Pm1 primer ( Figure 3A). In this figure, mir-6165 expression in Daoy, 1321N1, U87 (brain tumor-derived) cell lines are compared to a non-glioma (a human teratocarcinoma, NT2) cell line. mir-6165 and its host gene (NGFR) showed higher expression level in the brain tumor-derived cell lines, compared to the NT2 cells. This data supports the correct prediction of mir-6165 at the first place and also points to the origin of the cells (glioma) which express this miRNA at a detectable level.
Further supporting evidences were provided through the detection of endogenous precursor and mature form of mir-6165 in the brain tumor biopsies ( Figure 3B). mir-6165 endogenous expression was detected using real-time PCR on seven brain tumor tissue samples. Expression data was normalized using U48 endogenous control, and using one of the lowest degree meningioma samples as a reference. Along with the clear detection of endogenous mir-6165 in these tissue samples, a strong correlation was observed between mir-6165 and NGFR gene expression, which was further confirmed with Pearson's test (p = 0.0065). In all of the high grade tissue samples, the level of NGFR and mir-6165 were higher than the low grade samples ( Figure 3B).
In a neurotrophin-dependent manner, NGFR has been suggested as an important regulator of glioma invasion [15]. In that report, NGFR-positive cells within the glioma tumor samples are more migratory than the NGFR-negative glioma cells. Consistent with that report, we have seen upregulation of NGFR in higher grades of brain tumors.
A possible molecular mechanism in which mir-6165 co-operates with NGFR toward glioma invasion remains unknown; however, analysis of putative targets of mir-6165 by DAVID software [47] showed .50% of the predicted targets genes are those expressed in brain related tissues ( Figure S1) among them Pkd1 and DAGLA genes were down regulated in tested tumor samples (data not shown).

Cell Death Effect of pre-mir-6165 Overexpression in U87 Cell Line
In order to examine the effect of pre-mir-6165 overexpression on the rate of apoptosis induction in the transfected cells, flow  cytometry was performed using propidium iodide (PI) and annexin. This method has been reliably used to show the involvement of miR-16, let-7a and miR-34a in apoptosis [57]. Compared to the mock transfected control cells, no statistically significant apoptosis induction was detected 34 hours post transfection in Hela cells using PI staining. However, transfected U87 cells showed more than 10% elevation in sub-G1 cell population (Compare figure 4A with 5A). P value for such increased cell death rate was calculated by Repeated Measures ANOVA analysis and the results were highly significant (p = 0.0009). Annexin test shows the early apoptosis rate of the transfected cells and is more sensitive than PI test. This test showed ,40% of the transfected U87 cells were in early apoptosis (p = 0.003) ( Figure 4B). Overall, our results suggest a pro-apoptotic role for mir-6165 in the U87 cells, but not in the Hela cell line ( Figure 5). Lack of significant apoptotic effect of mir-6165 in transfected Hela cells may be due to the cellular environment or context of the cell [58] or lack of its partner(s) or target(s) in this cell line.
Accumulative experimental evidences from the detection of endogenous mature mir-6165 in the brain tumor samples to the detection of exogenous mir-6165 in the transfected Hela cells and its overexpression effect on the cell death rate, all emphasis the functionality of this miRNA.

Down Regulation of Predicted mir-6165 Target Genes upon its Overexpression
The sequence of mir-6165 was used as query for target prediction in DIANA-microT online tool and its top ten candidate genes are listed in Table 2. In the glioblastoma cell lines in which NGFR and mir-6165 are expressed, some of these target genes are down regulated based on the information from Oncomine database. DAGLA gene (Neural stem cell-derived dendrite regulator) is the highest scored target gene of the list and its miRNA recognition element (MRE) is highly conserved among other organisms. The interaction between this target and mir-6165 was further analyzed using RNAHybride (http://bibiserv.techfak. uni-bielefeld.de) online tool. Strong complementation was observed between the mir-6165 and its target MRE element in the DAGLA gene.
For experimental verification of the mir-6165 target genes, realtime specific primers were designed for (DAGLA and PKD) predicted target genes. Overexpression of pre-mir-6165 in Hela cells ended in 200,000 fold increase of mature form, 34h post transfection but DAGLA and PKD expression level still remained undetectable. In contrast to Hela cells, DAGLA and PKD genes were expressed at a detectable level in the untransfected U87 cells. Upon overexpression of mir-6165 in the U87 cells, a significant down regulation of these target genes was observed ( Figure 6).
Mutations in the PKD1 gene (encoding Polycistin-1) are accounted for the renal cysts formation [59]. Overexpression of PKD1 in Madin-Darby canine kidney (MDCK) cells has led to decreased apoptosis [60] and silencing of PKD1 has led to an increased apoptosis rate due to a reduced cell adhesion [61]. PKD1 and PKD2 are expressed in a number of tissues and organs, including the ductal epithelial cells in the kidney, liver, pancreas, breast, smooth muscle, endothelial cells of the vasculature and astrocytes in the brain [62]. U87 cell line is an astrocytoma cell line derived from a human malignant glioma and we showed that PKD1 was expressed in this cell line ( Figure 6). Consistent with previous studies as well as considering high expression level of PKD1 in most of the brain tumor-derived cell lines, a reduction in PKD1 expression following the overexpression of mir-6165 could justify the increased rate of apoptosis and a change of cell number distribution toward sub-G1 (Figure 4).
DAGLA gene is involved in endocannabinoid pathway shown by mouse knock out model. This pathway refers to a group of neuromodulatory lipids and their receptors [63]. DAGLA gene effect on proliferation is not yet clear to our knowledge. A recent publication has suggested that DAGLA deletion ended in ,50% less cell proliferation of the hippocampus cells [63]. DAGLA gene has multiple target sites for mir-6165 within its 39-UTR and it is deducible that overexpression of this miRNA has caused down regulation of this target gene at least at the mRNA level ( Figure 6).
It is remained to be tested if increased rate of apoptosis is related to the down regulation of DAGLA gene as well.
Co-expression of mRNAs and/or functional similarity of target and host genes of an intronic miRNA have been hypothesized [64]. This hypothesis was tested for NGFR and mir-6165 predicted target genes using GENEMANIA database ( Figure  S2). There was a co-expression pattern between most of the targets and NGFR host gene while none of the high scored targets show a direct physical interaction with NGFR gene.
It has been suggested that intronic miRNAs tend to target the genes that are functionally similar to their host genes [64]. Using Funsimmat algorithm (funsimmat.bioinf.mpi-inf.mpg.de) high scored target genes of mir-6165 showed functional similarity with NGFR gene as well. Noteworthy, Diana-mirpath and geneset2-miRNA online servers showed that most of the first 20 high scored targets of mir-6165 are also targeted by other miRNAs involved in cancers and brain development like mir-608, mir24 and mir-637. How this network of miRNAs with overlapped functions act in cancerous cell and development remain to be tested.
In conclusion, by multiple experimental evidence, our study revealed that mir-6165 is expressed in brain tumor-derived cell lines and primary brain tumor tissues. Overexpression of mir-6165 in U87 cell line increased the cell death rate and down regulated PKD1 and DAGLA target genes, which are involved in apoptosis. Furthermore, there is a strong co-expression network of mir-6165 host and target genes. Figure S1 Co-expression network of NGFR gene with mir-6165predicted targets. Although there is a co expression network between targets genes of mir-6165 but none of these targets have direct interaction with NGFR host gene. The targets of mir-6165have been shown by the larger nodes. (TIF) Figure S2 Percent of putative targets of Hsa-mir-6165 in different tissues. Analysis of mir-6165 putative target genes by DAVID showed more than 50 percent of predicted target genes are expressed in brain related tissues. (TIF)