PIWI-Interacting RNAs in Gliomagenesis: Evidence from Post-GWAS and Functional Analyses

An emerging area of cancer research concerns the putative involvement of PIWI-interacting RNAs (piRNAs), which are small (mostly 26–32 nt) noncoding RNAs with millions of members identified in mammalian cells (1). In humans, more than 30,000 piRNAs have been found, some of which map to hundreds or thousands of retrotransposable element–encoding genomic loci (2–4). piRNAs are best characterized by their role in guiding PIWI proteins and associated chromatin-silencing machinery to transposon-encoding DNA sequences in the genome (5–7). This mechanism results in DNA methylation of piRNA targets in mammals, thereby repressing the activity of potentially destabilizing transposable elements in a heritable manner (5). We and others have suggested that a subset of piRNAs may also be capable of regulating protein-coding genes via DNA methylation, which may bear on cancer development if the regulatory targets are cancer-relevant (8–10).

Recent evidence suggests that the PIWI–piRNA system may also be active in certain somatic tissues and that its dysregulation or reactivation may be a factor in carcinogenesis (11, 12). PIWI family proteins have been shown to be aberrantly expressed and associated with poor prognosis in ovarian, colorectal, gastric, hepatocellular, esophageal, and pancreatic cancers (13–18), and piRNA expression has been observed in bladder, colon, lung, cervical, breast, hepatocellular, and gastric cancers (19–24). Our recent study also identified a genetic variant in a mature piRNA sequence that is associated with breast cancer risk (25). In the case of glioma, it has been shown that PIWI family gene PIWIL1 displays oncogenic features in in vitro and in vivo glioma models (26) and that expression of the gene is associated with poor prognosis (27).

However, no study has explored the role that piRNAs may play in glioma development, despite the observation that some piRNAs are expressed in normal mammalian brain specimens (28). The objective of this study was therefore to examine whether inherited variants in piRNA-encoding loci are associated with glioma risk using data from 1,840 glioma cases and 2,401 controls from the GliomaScan Cohort (29) and to examine the underlying basis for one of the observed associations using in vitro functional studies.

Individual-level genotype data and phenotypic subject characteristics for participants of the GliomaScan Cohort–based Genome-wide Association Study (29) were downloaded from the Database of Genotypes and Phenotypes (dbGaP, Study Accession phs000652.v1.p1) after receiving data access authorization. Subjects were restricted to those genotyped on the Illumina Human660W-Quad BeadChip to avoid potential bias introduced by the case–control imbalance of subjects genotyped on other platforms. After removal of related or duplicate subjects (n = 25 cases, n = 150 controls), there were 1,840 cases (ICDO-3 codes 9380-9480) and 2,401 controls included for final analysis.

Genomic coordinates for all experimentally observed human piRNAs were obtained from the piRNABank database (30). SNPs included in the 1,000 Genomes Project Phase 3 (31) reference variant set (n = 77,818,332 biallelic SNPs) were identified within these coordinates. SNPs in piRNAs that map to more than 1,00 genomic loci were excluded, as these piRNAs may be less likely to be involved in protein-coding gene regulation (8).

Genotype and phenotype data were downloaded to a secure server at Yale University and decrypted and extracted according to dbGaP guidelines. 1,000 Genomes Phase III haplotypes were downloaded for use as the reference panel for imputation using IMPUTE v2.3.1 software (32). Input data were restricted to SNPs with call rate ≥ 90% and Hardy–Weinberg equilibrium HWE P > 0.0001 using the PLINK toolset (33). Fine mapping was conducted via imputation of all SNPs with minor allelic frequency (MAF) > 1% in 5-Mb chunks, and regional annotations were derived from the UCSC Genome Browser (34).

ORs and 95% confidence intervals (CI) for variant-glioma associations were estimated with SNPTEST v2.5 software (35), applying an additive allelic logistic regression model adjusting for sex, age, study design, and the first two principal components. Principal components were generated by the smartPCA algorithm in the EIGENSOFT v6.0 population genetics package to adjust for residual population substructure (36). Associations surpassing a Bonferroni-corrected significance threshold were deemed statistically significant and associations yielding false discovery rate-adjusted P < 0.10 were considered to be modest associations of interest (37).

Glioma cell lines U87 and A172 were purchased from ATCC and immortalized normal human astrocytes (NHA) were purchased from the University of California, San Francisco Tissue Core. All cell lines are tested for contaminants and authenticated prior to shipment using several characterization approaches including morphology assessment, karyotyping, cytochrome c oxidase variant analysis, and short tandem repeat DNA fingerprinting. Cells were not reauthenticated as they were passaged in our laboratory for fewer than 6 months after resuscitation. Cells were maintained in Eagle's modified essential medium (U87) or MEM (A172, NHA) and supplemented with 10% FBS. Single-stranded piRNA mimics were purchased from IDT (Supplementary Table S1), and single-stranded nontargeting RNA sequences of similar size were purchased from QIAGEN to act as negative controls in in vitro experiments. Cells were reverse-transfected using LipofectAMINE RNAiMAX transfection reagent (Invitrogen) according to the manufacturer's instructions.

Total RNA was isolated from U87, A172, and NHA cell lysates using the miRNeasy Mini Kit (QIAGEN), and cDNA was generated using the NCode miRNA First Strand cDNA Synthesis Kit (Invitrogen). qPCR was performed on an ABI-7500 System (Applied Biosystems) using a SYBR FAST qPCR Kit (Kapa Biosystems). Amplification reactions were conducted in triplicate with custom short piRNA forward primers and a universal reverse primer targeting appended poly(A) tails (Supplementary Table S1). Expression levels were normalized to small nuclear RNA U6 expression. Predicted secondary structures of piRNAs were generated by the Mfold v3.6 RNA folding algorithm using default parameters (38).

U87 cells were reverse transfected with either wild-type piR-598 or a nontarget control (NC) RNA. Cells were harvested 24 hours after transfection, and total RNA was isolated and approximately 1 μg was submitted to the Yale Center for Genome Analysis for genome-wide expression profiling on the Illumina HumanHT-12 v4 Expression BeadChip platform in biologic duplicate. Genes showing expression level differences between NC- and piR-598-WT treatments beyond a significance threshold of false discovery rate (FDR)-adjusted P = 0.05 were considered to be differentially expressed. Five genes were selected for expression validation by qPCR with input normalization to GAPDH. Network and pathway analyses were conducted using Ingenuity Pathway Analysis (IPA) software; P values for affected functional pathways were calculated using a Fisher exact test for enrichment of affected genes with a particular functional annotation. Array data have been deposited in the NCBI Gene Expression Omnibus repository (accession number GSE78935).

Cell viability was evaluated in piRNA-598- and negative control RNA-treated cell populations using the CellTiter 96 AQueous One Solution Cell Proliferation Assay (MTS) Kit (Promega). Briefly, cell viability was quantified at 48 and 96 hours after transfection. Color development was evaluated 1 hour after addition of MTS using a microplate spectrophotometer at an absorbance of 490 nm. Viability differences were analyzed using a Student's t test using six replicates per condition.

U87 cells were reverse-transfected with piRNA-598 or negative control oligos. Twenty-four hours after transfection, cells were trypsinized and resuspended in warmed EMEM with 0.36% agar. The mixture was plated in 60-mm cell culture dishes above a presolidified base layer of 0.75% agar. Dishes were incubated at 37°C with the addition of 500 μL complete media every 5 days. After 3 weeks, colonies were stained with 0.04% crystal violet/2% ethanol in PBS and photographed. Colonies were counted using ImageJ v1.48 software and compared between conditions using a Student's t test. Experiments were performed in triplicate.

To determine whether inherited variants in piRNA-encoding sequences are associated with the risk of adult glioma development, a genetic association analysis was conducted in 1,840 subjects with newly diagnosed glioma and 2,401 cancer-free controls included in the GliomaScan consortium. Approximately 67% of cases were diagnosed with high-grade glioma (grades III or IV), the majority of whom (82%) were of the glioblastoma multiforme subtype, and 55.2% of subjects were male (Supplementary Table S2).

A flowchart illustrating the selection of piRNA-sequence SNPs included in this analysis is presented in Supplementary Fig. S1. Of 2,514 SNPs of interest in piRNA-encoding sequences, 31 (1.2%) were directly genotyped on the Illumina HumanHap660W platform, and genotypes at 1,397 (55.6%) were imputed; 1,086 SNPs (43.2%) were excluded because they were unable to be imputed with sufficiently high quality due to low array coverage in piRNA-encoding intergenic regions. In total, 1,428 SNPs were analyzed for association with glioma risk adjusting for sex, age, study design, and the first two principal components. No evidence of systematic bias from underlying population substructure or other factors was detected in the input genotype data using this model (genomic inflation factor, λ = 1.009; Supplementary Fig. S2).

A Manhattan plot illustrates all 1,428 piRNA SNP–glioma association test results (Fig. 1A). Our analysis revealed a Bonferroni-corrected (P < 0.05/1,428 SNPs = 3.50 × 10⁻⁵) statistically significant association between glioma risk and rare variant rs149336947 (P = 2.34 × 10⁻⁵; FDR-adjusted P = 0.033), located near the 3′ end of piR-2799 on chromosome 2q33.1. piR-2799 is a 30-nucleotide piRNA that maps to the fourth intron of apoptosis inhibitor CFLAR, which is widely expressed in the human body including in the brain (39).

Four additional modest associations of interest were observed at rs62435800 in piR-18913 on chromosome 6q27 (P = 1.13 × 10⁻⁴; FDR-adjusted P = 0.054), rs147061479 in piR-598 on chromosome 8q13.1 (P = 1.69 × 10⁻⁴; FDR-adjusted P = 0.060), rs142742690 in piR-11714 on chromosome 9q22.1 (P = 1.10 × 10⁻⁴; FDR-adjusted P = 0.079), and rs35712968 in piR-3266 on chromosome 10q24.2 (P = 3.11 × 10⁻⁴; FDR-adjusted P = 0.089; Table 1).

To examine these associations at higher resolution, genotypes were imputed for all SNPs in 300-kb regions surrounding the five associated SNPs. For piR-2799, this analysis revealed nearly 100 SNPs with associations of comparable magnitude to that of rs149336947 spanning a 250-kb region of linkage disequilibrium (Fig. 1B). This region contains four genes and is upstream of one gene. In contrast, clusters of SNPs showing enhanced association signals were observed in more narrow regions of linkage disequilibrium surrounding rs62435800 in piR-18913, rs147061479 in piR-598, rs142742690 in piR-11714, and rs35712968 in piR-3266 (Fig. 1C–F). Both piR-18913 and piR-598 map to genetic regions that encode a small number of piRNAs and are devoid of protein-coding genes. piR-11714 is located on a piRNA-dense haplotype that does not encode protein-coding genes and is approximately 50 kb upstream of SPIN1, and piR-3266 maps to the 3′ untranslated region (UTR) of HPSE2.

piRNA expression measurement was conducted using a qPCR-based method. Results showed expression of piR-18913, piR-598, piR-11714, and piR-3266 in all cell lines tested (U87, A172, and NHA; Supplementary Fig. S3). Expression of piR-2799 was not detected in any of the cell lines.

Among the four candidate piRNAs (piR-18913, piR-598, piR-11714, and piR-3266) that were found to be expressed in the cell lines examined, piR-598 harbored the variant conferring the greatest magnitude of glioma risk or protection and therefore was the subject of additional in vitro functional analyses. The predicted secondary structure of piR-598 is illustrated in Fig. 2A. In the most thermodynamically stable structure, the piRNA forms a small hairpin loop from the 5th to 19th bases; variant rs147061479, located at the 29th of 31 bases, is not involved in the predicted loop structure. We performed transcriptome-wide expression profiling 24 hours after transient upregulation of the piRNA in U87 cells to examine the impact of the expression of this piRNA in the context of glioma. Relative to nontargeting control-treated cells, a total of 518 transcripts were observed to be differentially expressed at FDR-adjusted P < 0.05 (Supplementary Table S3), the majority of which (71.2%) were observed to be underexpressed in piR-598–treated cells. Expression differences for five transcripts selected for validation of expression array data by qPCR were generally consistent with array results (Supplementary Fig. S4).

Subsequent pathway analysis showed that piR-598–affected genes were significantly enriched for those involved in cell death and survival (P = 3.43 × 10⁻³), cell-cycle progression (P = 2.63 × 10⁻³), and cellular assembly and organization (P = 2.39 × 10⁻³; Fig. 2B). Network visualization analysis revealed a core of functionally interrelated molecules including BAX, a key regulator of p53-mediated apoptosis, and oncogenic transcription factor JUN (Fig. 2C).

Wild-type and variant piR-598 mimics were independently overexpressed in glioma (U87 and A172) and NHA cell lines, and cell viability was measured. Relative to a nontargeting control RNA, transfection of wild-type piR-598 sharply reduced proliferation of both U87 and A172 cells, notably with nearly 40% inhibition measured 96 hours after transfection in U87. However, transfection of the mutant rather than wild-type piR-598, containing the variant allele, significantly attenuated the antiproliferative impact. The same pattern was observed in normal glial cell line NHA (Fig. 3A).

The functional impact of the piRNA variant was also examined on U87 colony formation in soft agar, which is a model of anchorage-independent growth potential. Treatment with wild-type piR-598 reduced the number of colonies formed to approximately half those formed following negative control treatment. However, treatment with the variant rather than wild-type piRNA was sufficient to not only eliminate the anti-proliferative effect of the piRNA but also confer a more than 4-fold increased colony-forming potential relative to wild-type piR-598 treatment (Fig. 3B).

This post-GWAS study suggests that inherited variants at five piRNA loci (FDR-adjusted P < 0.10) are associated with glioma risk in the GliomaScan Cohort, the largest publicly available glioma GWAS dataset. None of these associations has been reported in previous publications. Genomic loci at 8q24 and 9p21 have been linked to glioma in previous GWAS (40, 41); however, observed associations on these chromosomes at rs147061479 (piRNA-598 at 8q13) and rs142742690 (piRNA-11714 at 9q22) are unrelated to previous signals, suggesting that novel genetic risk loci harboring piRNA variants have been identified from our post-GWAS approach.

A Bonferroni-adjusted significant association was detected between glioma risk and rs149336947 in piR-2799. Regional imputation showed that this association extended over a large region of linkage disequilibrium that harbors four genes including apoptosis regulator CFLAR, as well as the promoter region for initiator caspases CASP10 and CASP8 that has been linked to susceptibility to several cancers (42). Thus, the association at rs149336947 may reflect a functional polymorphism that is unrelated to piR-2799, possibly representing a separate low-frequency biomarker of glioma risk that is itself worthy of further follow-up. The observation that piR-2799 expression was undetectable in U87, A172, or NHA cell lines further supports this notion.

In contrast, regional imputation of the other four regions harboring piRNAs piR-18913, piR-598, piR-11714, and piR-3266 revealed narrow clusters of SNPs with amplified association signals relative to surrounding areas. There are no protein-coding genes in the regions that encode piR-598 and piR-18913. Moreover, expression of all four of these piRNAs was detectable in both normal glial and glioma cell lines. These findings suggest potential biologic roles of these piRNAs and their variants in gliomagenesis that warrants further examination.

The functional significance of one of the identified piRNAs, piR-598, was explored in follow-up transcriptional profiling and network analyses, which indicated involvement of this piRNA in cell death and survival pathways. Of particular interest was the upregulated expression of the BAX transcript, which encodes a protein that promotes cell death by inhibiting apoptosis repressor Bcl-2 (43). Expression of the closely related GOS2 gene, encoding another Bcl-2–interacting and apoptosis-promoting protein (44), was upregulated. We performed expression profiling after a relatively short treatment period (24 hours) to detect early piRNA-induced transcriptional changes before cell viability was compromised; gene expression differences in this experiment did not tend to be large in magnitude as a result.

Subsequent in vitro assays confirmed the role of piR-598 in cellular growth identified from the expression profiling analysis and further demonstrated the functional impact of the genetic variant. Delivery of the wild-type piRNA-598 mimic significantly diminished cell viability relative to control treatment. However, upregulation of the variant rather than wild-type piR-598 sharply attenuated the antiproliferative response observed with the wild-type piRNA. Additional evidence comes from the observation that wild-type piR-598 treatment limited long-term colony formation of U87 cells seeded in soft agar, yet treatment instead with the variant piRNA was sufficient to eliminate the antiproliferative effect of piR-598 and in fact promoted colony formation. The discrepancy in the effect of the variant piRNA with respect to negative control treatment in the two assays was likely attributable to the difference in time period of cell growth after piRNA treatment (4 vs. 21 days), as an increased growth rate attributable to the variant piRNA was more readily revealed via the greater number of population doublings occurring in the longer term colony formation assay. These results provide consistent functional support for the increased glioma risk associated with rs147061479.

While these post-GWAS and functional analyses suggest a role of PIWI-interacting RNAs in gliomagenesis, the current study has some limitations. Although the association at rs147061479 in piR-598 is supported by follow-up functional analyses, we recognize that replication of the piR-598 association in an independent population would provide additional support for its role in gliomagenesis. Furthermore, because of the unavailability of individual-level glioma subtype data, we were unable to explore subtype-specific associations despite the widely accepted heterogeneity of the disease. However, in vitro analyses of piR-598 were performed in cell line models of grade IV glioma, supporting a role for the rs147061479 variant in this subtype. In addition, genotypes of glioma-associated piRNA SNPs were imputed rather than directly genotyped; however, we utilized a strict threshold for imputation quality and accounted for genotype uncertainty in our model, thus resulting in appropriately conservative association estimates. Finally, additional mechanistic work will be required to uncover the precise molecular basis for the observed genetic associations and phenotypic impact of piR-598. We and others have demonstrated that piRNAs may be capable of inducing DNA methylation of protein-coding genes (8, 10), and PIWI family protein PIWIL4 has been shown to localize to the nucleus and may be involved in de novo DNA methylation of transposable elements (45) as well as targeted protein-coding genes. It has also been suggested that certain PIWI–piRNA complexes may have siRNA-like roles as posttranscriptional regulators in mRNA degradation and translational regulation (1, 46). Establishing the direct targets of piRNAs in a cancer context is an important area for future investigation.

In summary, this multidisciplinary study demonstrates a novel role of piRNAs, a group of small noncoding RNAs, in gliomagenesis using evidence from both post-GWAS and in vitro functional analyses. These findings add to a growing body of research linking the PIWI–piRNA pathway and human cancer and strongly support expanded exploration of this abundant yet relatively understudied class of noncoding RNA in relation to cancer development. Pursuing this research may yield novel biomarkers for cancer risk, insights into the tumor developmental process, and strategies for cancer treatment.

No potential conflicts of interest were disclosed.

Conception and design: D.I. Jacobs, A. Fu, R. Dubrow, H. Lin, Y. Zhu

Development of methodology: D.I. Jacobs, Q. Qin, Y. Zhu

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D.I. Jacobs, Q. Qin, G. Wang, Y. Zhu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D.I. Jacobs, A.T. DeWan, H. Lin, Y. Zhu

Writing, review, and/or revision of the manuscript: D.I. Jacobs, Q. Qin, A. Fu, R. Dubrow, E.B. Claus, A.T. DeWan, G. Wang, H. Lin, Y. Zhu

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Zhu

Study supervision: E.B. Claus, H. Lin, Y. Zhu

This work was supported by the NCI grant CA122676 and funds from Yale University. D.I. Jacobs was supported by the Yale–NCI predoctoral training grant T32 CA105666. A. Fu was partially supported by UCLA-NCI training grant T32 CA09142. E.B. Claus was partially supported by CA119215.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.