Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Schizophrenia-Associated MIR204 Regulates Noncoding RNAs and Affects Neurotransmitter and Ion Channel Gene Sets

  • Sophia Cammaerts,

    Affiliations University of Antwerp, Antwerp, Belgium, Applied Molecular Genomics Unit, Department of Molecular Genetics, VIB, Antwerp, Belgium

  • Mojca Strazisar,

    Affiliations University of Antwerp, Antwerp, Belgium, Applied Molecular Genomics Unit, Department of Molecular Genetics, VIB, Antwerp, Belgium

  • Bart Smets,

    Affiliations University of Antwerp, Antwerp, Belgium, Centralized Service Facility, Department of Molecular Genetics, VIB, Antwerp, Belgium

  • Sarah Weckhuysen,

    Affiliations University of Antwerp, Antwerp, Belgium, Neurogenetics Group, Department of Molecular Genetics, VIB, Antwerp, Belgium

  • Annelie Nordin,

    Affiliation Division of Psychiatry, Department of Clinical Sciences, Umeå University, Umeå, Sweden

  • Peter De Jonghe,

    Affiliations University of Antwerp, Antwerp, Belgium, Neurogenetics Group, Department of Molecular Genetics, VIB, Antwerp, Belgium, Antwerp University Hospital, Antwerp, Belgium

  • Rolf Adolfsson,

    Affiliation Division of Psychiatry, Department of Clinical Sciences, Umeå University, Umeå, Sweden

  • Peter De Rijk,

    Affiliations University of Antwerp, Antwerp, Belgium, Applied Molecular Genomics Unit, Department of Molecular Genetics, VIB, Antwerp, Belgium

  • Jurgen Del Favero

    jurgen.delfavero@molgen.vib-ua.be

    Affiliations University of Antwerp, Antwerp, Belgium, Applied Molecular Genomics Unit, Department of Molecular Genetics, VIB, Antwerp, Belgium, Multiplicom N.V., Niel, Belgium

Abstract

As regulators of gene expression, microRNAs (miRNAs) are likely to play an important role in the development of disease. In this study we present a large-scale strategy to identify miRNAs with a role in the regulation of neuronal processes. Thereby we found variant rs7861254 located near the MIR204 gene to be significantly associated with schizophrenia. This variant resulted in reduced expression of miR-204 in neuronal-like SH-SY5Y cells. Analysis of the consequences of the altered miR-204 expression on the transcriptome of these cells uncovered a new mode of action for miR-204, being the regulation of noncoding RNAs (ncRNAs), including several miRNAs, such as MIR296. Furthermore, pathway analysis showed downstream effects of miR-204 on neurotransmitter and ion channel related gene sets, potentially mediated by miRNAs regulated through miR-204.

Introduction

miRNAs are short regulatory noncoding RNA molecules. Their maturation process consists of two consecutive cleavage steps. The first step takes place in the nucleus, whereby the Microprocessor complex, consisting of ribonuclease Drosha and its cofactor DGCR8, frees the hairpin structured precursor miRNA (pre-miRNA) from the primary miRNA transcript (pri-miRNA) [1,2]. After export to the cytosol by Exportin-5 [3], the ribonuclease Dicer cuts the terminal loop off the pre-miRNA [4,5]. Subsequently, the mature miRNA duplex, containing the 5p and 3p mature miRNA strands, is loaded onto the RNA-induced silencing complex (RISC). One of the strands is retained and will guide Argonaute, the active compound of RISC, to its targets via incomplete sequence complementarity [6,7]. Posttranscriptional regulation of messenger RNAs (mRNA) in the cytosol by miRNAs is well studied [8]. Next to this canonical cytoplasmic role of miRNAs, mature miRNA-Argonaute complexes can be shuttled back into the nucleus by Importin-8 [9]. Many mature miRNAs have been shown to be present in the nuclear as well as the cytoplasmic fraction of cells [10], indicating functionality in the nucleus. Recently, it was also shown that though mRNAs are the main targets of miRNA regulation, as much as 30% of miRNA interactions take place with other classes of RNAs [11].

Identification of miRNAs involved in disease pathomechanisms is often performed at the RNA level, by comparing miRNA expression profiles in affected tissues between individuals with different phenotypes or in affected versus nonaffected tissue of patients [12,13]. The miRNAs that are differentially expressed between the studied states then indicate which miRNAs may be functionally relevant for the disease. However, from such analysis it is not possible to derive whether the miRNA is a primary factor in the process or a downstream effect [14]. In addition, for diseases where the affected patient tissue is not readily available, the profiling approach is not applicable. Another starting point is a genetic approach, where one searches for causal or risk factors associated with disease in (miRNA) genes. The miRNA genes with associated or causal variants can then be further investigated to assess the functional effects of these variants [15,16].

Here we describe a high throughput genetic approach for the identification of miRNAs regulating tissue-specific functions. This strategy is based on a large-scale targeted miRNA gene variant screening and association analysis with disease phenotypes related to the tissue of interest. Variants enriched in these phenotypes are then functionally investigated to study the wild type function of these miRNA genes. We applied this strategy to identify miRNAs with a function in neuronal pathways by variant discovery in brain expressed miRNA genes in patients with neurological disorders: schizophrenia and idiopathic generalized epilepsy (IGE). Two selected miRNA gene variants were functionally validated in the neuronal-like cell line SH-SY5Y.

Results

Large-scale variant screening

To study the function of miRNA genes in neurological diseases and to investigate which miRNA genes regulate neuronal processes, we created a Multiplex Amplification of Specific Targets for Resequencing (MASTR) assay [17,18] for the amplification of 289 brain expressed miRNA genes. This assay was subsequently used in combination with massively parallel sequencing for variant discovery within the selected miRNA genes in patients with schizophrenia or IGE and in Swedish and Belgian and Dutch control individuals. After mapping, variant calling and variant filtering, an association analysis was performed to assess which variants are over- or underrepresented in the disease phenotypes and thus may have a neuronal function (S1 File). After filtering, a total of 265 variants were identified in the schizophrenia patients and the Swedish control group and 315 variants were identified in the IGE patients and the Belgian/Dutch control group (Tables A and B in S1 File). In both patient groups, individual patients had on average 45 different variants in or near miRNA genes (Figure A in S1 File). In the schizophrenia group, ten variants had a difference in variant frequency between the patients and controls (unadjusted p < 0.05). One of these variants, rs7861254, was significantly associated with schizophrenia after Bonferroni correction. In the IGE group, eight variants had a difference in frequency in the patients compared to the controls (Tables C and D in S1 File). All variants with (unadjusted) p-value < 0.05 were considered for assessment of variant location relative to the functional miRNA regions and structural impact prediction using the miRVaS software [19].

Based on the association analysis and impact prediction, two variants were selected for further functional studies. The first variant, rs7861254 C>T, showed a significant lower variant frequency in schizophrenia patients compared to Swedish control individuals (alternative allele frequency patients / controls: 16.39% / 25.34%, unadjusted p = 1.51x10-4; Bonferroni adjusted p = 0.04). Although the variant is located 107 nt outside the MIR204 hairpin, it is predicted to induce structural changes in the hairpin end and flanking regions, which can have an effect on its processing (Fig 1A and 1B). The second variant, rs2682818 A>C, is more frequently present in IGE patients than in Belgian/Dutch controls (alternative allele frequency patients / controls: 91.67% / 85.75%, unadjusted p = 9.63x10-3, Bonferroni adjusted p = 3.03) and is located in MIR618, in the pre-miR-618 hairpin arm opposite of the mature miR-618 sequence. Based on the secondary structure prediction, we hypothesized that the variant allele may affect the distribution pattern of canonical and alternative mature miR-618 sequences (isomiRs), as it was predicted to induce a structural change at the site of Drosha cleavage (Fig 1C and 1D). This structural change was very consistent between different types of structure predictions (minimal free energy (MFE), centroid, maximal expected accuracy (MEA)) and with varying flank sizes.

thumbnail
Fig 1. MFE secondary structure predictions of hsa-mir-204 and hsa-mir-618 as generated by miRVaS.

(a-b) MFE structure for hsa-mir-204 with variant rs7861254, (a) wild type miRNA, (b) variant miRNA. As the variant is located 107 nt outside the hairpin, flanking regions of 150 nt were included for the predictions. Centroid and MEA predictions with this flank size also showed large structural changes (but different changes), while predictions with flanks of 200 nt resulted in minor changes far away from the hairpin. (c-d) MFE structure of hsa-mir-618 for variant rs2682818, (c) wild type miRNA, (d) variant miRNA. The variant is predicted to induce a shift of the first base of the miR-618 sequence into an internal loop. Flanks of 100 nt were used for this prediction. Predictions were also run for the hairpin with flank sizes of 50 nt, 150 nt and 200 nt and centroid and MEA structures: all predicted the same change within the hairpin. Color scheme: magenta: mature miRNA, orange: seed region, dark blue: terminal loop, cyan: hairpin. The variant is colored in red and indicated by an arrow. Structural changes induced by the variant are colored in black.

https://doi.org/10.1371/journal.pone.0144428.g001

Functional characterization of MIR204

Variant rs7861254 results in altered processing of miR-204.

To explore the functional effect of rs7861254 genotypes and to investigate MIR204 function, we established stable cells overexpressing the wild type or mutant MIR204 gene. As our aim was to study this gene in the context of neurological disease, we used the neuronal-like cell line SH-SY5Y for these experiments. We assessed the expression of both mature miRNAs (miR-204-5p and miR-204-3p) in replicates of wild type (MIR204WT) and mutant (MIR204SNP) cell lines by RT-qPCR and corrected expression for transduction efficiency. In a first set of triplicates of wild type and mutant cells, a lower expression of both mature miRNA sequences was measured in the mutant cells (miR-204-5p: fold change MIR204SNP versus MIR204WT = 0.166, p = 0.046; miR-204-3p: fold change MIR204SNP versus MIR204WT = 0.165, p = 0.050; Fig 2A, Figure B(a) in S1 File), confirming the predicted effect of the variant on the biogenesis of miR-204. A comparable trend was found in a set of three additional clones (Figures B(b) and C in S1 File, p > 0.05).

thumbnail
Fig 2. QPCR (a), microarray (b) and small RNA sequencing results (c) for MIR204 cells.

(a) Average normalized relative quantity (NRQ) for miR-204-5p and miR-204-3p in MIR204WT (CC genotype, black bars, clones 1–3) and MIR204SNP cells (TT genotype, white bars, clones 1–3) normalized for transduction efficiency. (b) Average normalized expression values (NEV) for significantly differently expressed genes between MIR204WT (CC genotype, black bars, clones 4–5) and MIR204SNP (TT genotype, white bars, clones 4–6) cells. (c) Log2 fold changes of miRNAs significantly differently expressed in MIR204SNP cells (clones 4–6) compared to MIR204WT cells (clones 4–5). Error bars in panels (a-b) represent standard deviation of biological replicates, error bars in (c) represent standard error as calculated by DESeq2.

https://doi.org/10.1371/journal.pone.0144428.g002

Transcriptome analysis.

We next analyzed the effect of the variant-induced lowered expression of miR-204 on the transcriptome of the stable SH-SY5Y cells using microarray analysis to deduce the wild type function of MIR204. Significance Analysis of Microarray (SAM) analysis indicated that 15 genes are significantly differentially expressed in the mutant cells versus the wild type cells (Fig 2B, Table 1, Figure D in S1 File). Rather unexpectedly, these genes all code for ncRNAs. The gene with highest fold change is MIR296, which is upregulated 33-fold in the MIR204SNP cells.

thumbnail
Table 1. Genes significantly differentially expressed between MIR204WT and MIR204SNP cells, as identified by SAM.

https://doi.org/10.1371/journal.pone.0144428.t001

To investigate the effect of miR-204 on whole pathways as opposed to separate genes, we also performed gene set enrichment analysis (GSEA). GSEA was performed with Gene Ontology gene sets for molecular function, comparing MIR204WT to MIR204SNP cells (Tables 2 and 3). 30 gene sets were significantly enriched in MIR204WT cells. The gene set with the strongest normalized enrichment score is neurotransmitter binding. In addition, several ion channel activity related gene sets are highly enriched. In MIR204SNP cells, 17 gene sets were significantly enriched. These sets were related to RNA and DNA binding and ligase activities, amongst other functions.

thumbnail
Table 2. Gene sets enriched in MIR204WT cells with FDR < 0.01 and nominal p-value < 0.001.

https://doi.org/10.1371/journal.pone.0144428.t002

thumbnail
Table 3. Gene sets enriched in MIR204SNP cells with FDR < 0.01 and nominal p-value < 0.001.

https://doi.org/10.1371/journal.pone.0144428.t003

Since several ncRNAs are significantly upregulated in MIR204SNP cells, it is plausible to assume a more direct relationship between the RNA, DNA binding and ligase gene sets and miR-204 function and an indirect relationship with the neurotransmitter and ion channel gene sets, possibly via miR-296 or any of the other ncRNAs that are differentially expressed in the mutant cells.

Genome-wide mature miRNA expression analysis.

To evaluate the effect of miR-204 on mature miRNA expression, the cells were also subjected to small RNA sequencing. Differential miRNA expression analysis using DESeq2 [20] indicated seven mature miRNAs that are significantly differently expressed between MIR204WT and MIR204SNP cells, in addition to miR-204-5p and miR-204-3p (Benjamini-Hochberg corrected p-value < 0.1; Fig 2C). Consistent with the microarray results where MIR296 is strongly upregulated, miR-296-5p and miR-296-3p are both significantly upregulated (9- and 12-fold) in MIR204SNP cells. Furthermore, consistent with the presence of other ncRNAs being upregulated according to the array results, five other mature miRNAs are deregulated (three upregulated, two downregulated), indicating that this might be a common function of MIR204.

Target prediction.

To assess whether the differentially expressed transcripts contain binding sites for miR-204-5p or miR-204-3p regulation, we predicted the presence of seed binding sites using the RNAhybrid software [21]. Because the primary transcript length of miRNAs is not known in most cases, we included a flanking region of 200 nt surrounding the miRNA hairpin. Of the 18 transcripts tested, nine contain at least one seed match, including pri-miR-296 (S5 Table). Of these, pri-let-7i has the highest number of potential binding sites (after correction for transcript length), two for miR-204-5p and one for miR-204-3p. Though the presence of a seed match can further strengthen the finding that they are regulated by miR-204, the absence of a seed match not necessarily weakens it since it has been shown recently that miRNA-target interactions without a seed match are quite common [11].

Functional characterization of MIR618

To investigate the effect of the variant rs2682818 and the function of MIR618, stable SH-SY5Y cells overexpressing wild type or mutant MIR618 were established and subjected to the same analyses as MIR204 cells. The expression level of miR-618 in the MIR618WT cells was too low to determine normalized relative quantities (NRQ) of miR-618 between MIR618WT and SNP cells by RT-qPCR. We reasoned this could either be due to very low or no expression of miR-618 in MIR618WT cells and/or due to the expression of an isomiR that cannot be detected by the assay used for RT-qPCR.

SAM analysis in MIR618 cells failed to identify genes that were significantly differentially expressed in the mutant cells compared to the wild type cells (Figure E in S1 File). GSEA showed three gene sets to be significantly enriched in MIR618SNP cells (Table 4): helicase activity, taste receptor activity and ATPase activity. No pathways were enriched in MIR618WT cells with FDR < 0.05.

thumbnail
Table 4. Gene sets enriched in MIR618SNP cells with FDR < 0.01 and nominal p-value < 0.001.

https://doi.org/10.1371/journal.pone.0144428.t004

Differential miRNA gene expression analysis on the small RNA sequencing data indicated that miR-618 was significantly differentially expressed with a 10-fold change. Because the secondary structure predictions indicated a potentially altered isomiR distribution, we investigated whether this was indeed the case. For this analysis, we used the raw read counts of all reads mapping to the miR-618 locus and grouped them in two classes: reads with the same 5' start as the canonical miRNA sequence as deposited in miRBase [22] (termed “5’ canonical isomiRs”) and reads with a different 5' start position (termed “5’ change isomiRs”), as these might have a different functionality. Per class the reads were summed and normalized to the total read count on the miR-618 locus within the sample. The MIR618SNP cells, that have a higher expression of miR-618, show a significantly reduced frequency of 5' change isomiRs compared to MIR618WT cells (p = 0.040; Figure F in S1 File). Given that the miR-618 expression in MIR618WT cells is very low (20–75 reads in total), the increased relative presence of 5' isomiRs may mean that these reads are degradation products rather than functional isomiRs. For comparison, we did the same isomiR analysis for miR-204-5p and miR-204-3p in MIR204 cells (Figure F(b and c) in S1 File. The frequency of 5' change isomiRs of miR-204 was not significantly different between MIR204 genotypes (miR-204-5p: p = 0.146, miR-204-3p: p = 0.573), indicating a MIR618 genotype-dependent specific event.

Given that miR-618 expression is upregulated in the MIR618SNP cells, the enriched gene sets are likely indirect targets of miR-618 since their expression is positively correlated with miR-618 expression. As no genes were significantly differently expressed after a 10-fold upregulation of miR-618, gene sets were only enriched in the SNP cells and the expression of miR-618 in MIR618WT cells was very low, this may indicate that miR-618 has no direct targets in this specific cell type, and that ectopic overexpression results in dysregulation of downstream pathways.

Discussion

We present a novel large-scale approach for the identification of miRNA genes involved in neuronal functions. The identification of relevant miRNA genes is based on an association study of variants in brain expressed miRNA genes in patient groups with neurological phenotypes and was performed here on schizophrenia and IGE. The genetic association of a miRNA variant with one of the phenotypes is employed as an indication for a miRNA to have a potential role in the regulation of neuronal pathways. This allowed us at the same time to identify genes with a potential function in neuronal processes and to identify variants that may be significantly associated with the studied phenotypes, such as the variant rs7861254 that is associated with schizophrenia. Functional validation of the selected miRNA genes was also performed on a large scale, implementing both genome-wide mRNA and miRNA profiling, to assess all possible targets of these genes and to assess whether the variants had any effect on the isomiR distribution pattern of the miRNA. This approach led to interesting and different findings for MIR204 and for MIR618. The strategy can be adapted to identify miRNAs involved in a wide range of other organ- or tissue-specific functions, provided phenotypes are known, where the tissue or organ of interest is affected, and are available for genetic variant screening.

MIR204 is located within an intron of TRPM3, which encodes for a transient receptor potential cation channel that plays a role in detecting noxious heat and chemical stimuli in sensory neurons [23]. Several cancer-related studies have investigated the function of MIR204. miR-204-5p was shown to have a tumor-suppressor function through the negative regulation of cell migration and invasion [2427]. Its expression is reduced in different cancers, such as endometrial cancer, glioma, retinoblastoma and colorectal cancer [2427]. miR-204-5p was demonstrated to regulate expression of TRPM3 both directly and indirectly in renal cell carcinoma cell lines, as part of a network involved in regulation of autophagy [28]. miR-204 also plays an important role in the eye. A causal mutation in the miR-204-5p seed sequence was recently found in a family with inherited retinal dystrophy and coloboma [29]. In vivo studies in medaka fish implicated miR-204 in photoreceptor differentiation and function [29]. Previously, miR-204 had already been shown to play a crucial role in development of the lens and retina in medaka fish [30]. Transcriptome analysis in whole embryos also showed nervous system development, neurogenesis and axon guidance pathways to be differentially regulated when miR-204 expression was manipulated [31]. Its role in axonal outgrowth of retinal ganglion cells was further validated in vivo. A general mode of action of miR-204-5p appears to be the regulation of transcription factors. For example, miR-204-5p has been shown to regulate a network of transcription factors, MEIS2, FOXC1, PAX6, involved in eye development [30,32] and it targets the insulin transcription factor MAFA in pancreatic beta-cells [33]. In addition, one third of known miR-204-5p targets (9/27) listed in miRTarBase 4.5 [34] (only taking into account those interactions directly validated using reporter assays) are members of different transcription factor families, confirming this to be an important function of miR-204-5p.

The variant rs7861254 has not yet been functionally investigated for its effect on miR-204 expression. The alternative allele is significantly less frequent in schizophrenia patients than in the control individuals. By functionally validating the effects of rs7861254, we show that miR-204 regulates other miRNAs and ncRNAs. It exerts its largest effect on miR-296: both the pri/pre-miR-296 (as measured on the array) as miR-296-5p and -3p (small RNA sequencing data) are significantly upregulated in the MIR204SNP cells that have reduced miR-204 expression. Functional validation of miR-618 on the other hand did not identify any significantly differentially expressed target genes, when applying the same approaches. To the best of our knowledge, a direct regulatory role for miR-204 of other miRNAs has not been described yet, though it is consistent with the known function of miR-204 to target transcription factors as it characterizes miR-204 as a regulator of different classes of gene expression regulators in the cell.

Next to mediating posttranscriptional regulation of mRNAs, miRNAs can target promoters of protein-coding target genes resulting in transcriptional activation or repression [35]. It is therefore likely that they can also mediate the same type of transcriptional regulation of miRNA genes [36], either by targeting transcription factors that regulate the expression of miRNAs, or by targeting miRNA genes directly. In addition, both maturation and stability of miRNAs have been shown to be regulated by other miRNAs: mouse specific miR-709 regulates miR-15a-16-1 biogenesis by targeting the primary transcript and preventing its processing to pre-miRNA [37], while base pairing between mature miR-107 and let-7 reduces their stability bidirectionally [38]. That these miRNA-miRNA interactions are not isolated events was shown by a large-scale analysis of direct miRNA-target interactions with the CLASH technology in HEK293 cells where 3% of all interactions occurred between miRNAs [11].

We did not find any transcription factors or other protein-coding genes to be significantly deregulated upon miR-204 expression manipulation and did find potential miR-204 binding sites in several miRNAs, pointing to a more direct role of miR-204 in regulation of other ncRNAs. Nevertheless, the possibility of intermediary targets, for instance those targets with altered protein expression but unaltered mRNA expression, cannot be ruled out completely. The upregulation of MIR296 expression at different levels indicates that the regulation occurs upstream of the pre-miRNA to mature miRNA processing and may take place either at pri-miRNA processing or at the level of transcription. For the other mature miRNAs that have different expression in MIR204SNP cells, the pri/pre expression is not significantly affected. Both a biological and a technical reason may exist. It could be that there are indeed no significant alterations at the pri/pre-miRNA level, pointing toward regulation at a different step compared to MIR296, e.g. Dicer processing or miRNA destabilization. On the other hand, it could be that the changes at the pri/pre level are not large enough to be picked up by the technology used. Further experiments are required to determine the direct or indirect nature of these interactions and at which step they take place.

MIR618 is located within the first intron of the gene LIN7A. LIN7 proteins are present in pre- and post-synaptic protein complexes and are involved in regulating neurotransmitter release and transport of NMDA receptor vesicles [3942]. The function of MIR618 has not been explored in detail yet. Overexpression of miR-618 in an anaplastic thyroid cancer cell line resulted in inhibition of proliferation, invasion and migration of these cells [43]. RIP-Chip analysis of HeLa cells transfected with a miR-618 mimic and subsequent network analysis indicated a potential role for MIR618 in lymphomagenesis [44]. In the same study, the A allele (wild type allele) of rs2682818 was found to be associated with an increased risk for follicular lymphoma. In vitro validation in HeLa cells showed that the AA genotype results in decreased miR-618 expression compared to the alternative CC genotype.

In this study, we found that the variant rs2682818 has a higher frequency in patients with IGE compared to control individuals (though not statistically significant after correction for multiple testing). Validation in SH-SY5Y cells indicated that the homozygous variant genotype (CC) results in a significantly higher expression of miR-618 in SH-SY5Y cells compared to the AA genotype. In addition, the variant cells have a significantly lower frequency of 5’ change isomiRs than MIR618WT cells, which confirmed our structure prediction-based hypothesis that the variant may induce altered pre-miR-618 processing resulting in a changed isomiR distribution pattern. Nevertheless, the large expression change induced by the variant did not lead to expression changes of protein-coding or of other microRNA genes. The absence of deregulated miRNA genes shows that the miRNA-regulation mode of action is a finding specific to miR-204. On the other hand, the presence of many 5’ changed isomiRs and the lack of enriched gene sets with negative correlated expression despite a large change in expression of miR-618 between MIR618WT and SNP cells may indicate that miR-618 has no endogenous role in the studied cell type and that high ectopic expression of this miRNA results in side effects.

With this study we demonstrated the different effects of two selected genetic variants on miRNA functionality. The variant rs2682818 in MIR618 resulted in increased miR-618 expression and in a different distribution of miR-618 isomiRs. The variant rs7861254 near MIR204 is significantly associated with schizophrenia, identifying a new instance of a miRNA involved in schizophrenia etiology. By studying this variant, that reduced miR-204 expression, a new function of MIR204 became apparent: the regulation of other ncRNAs. Though the mechanism of miRNAs regulating expression of other miRNAs has been described before, it remains as yet an underexplored functionality of miRNAs. The discovery of a new schizophrenia related miRNA gene regulating neurotransmitter and ion channel gene sets and with a novel mode of action shows the utility of our large scale miRNA gene sequencing and functional analysis approach and indicates that it will be useful for the study of miRNAs in other tissues as well.

Methods

Ethics statement

The study was approved by the Medical Ethical Committees of the universities of Umeå and Antwerp. Each individual signed an informed consent form for participation in the study.

Samples

DNA of patients with schizophrenia, of patients with IGE and of control individuals was collected. The schizophrenia patients included in this study were recruited in Sweden (n = 186). Control individuals for this patient group were also recruited in Sweden (n = 1043), as part of the Betula study [45,46]. Patients with IGE were recruited in Belgium (n = 162). Control individuals for this patient group were recruited in Belgium and the Netherlands (n = 273).

Variant screening

A MASTR assay that amplifies 289 brain expressed miRNA genes was created (assay available upon request). This assay was subsequently used to screen the genomic DNA of patients and control individuals. Control DNA samples were pooled per seven samples prior to amplification, patient samples were amplified individually. Libraries were sequenced on MiSeq (Illumina, CA, USA) using v3 chemistry. After performing filtering steps in GenomeComb [47] (S1 File), remaining variants were tested for association with the schizophrenic or epileptic phenotype using a two-sided Fisher’s exact test. All variants with p < 0.05 were subsequently analyzed using the miRVaS software [19] to assess their location in the miRNA (seed, mature, terminal loop, hairpin, flanks) and their predicted structural impact. Bonferroni adjusted p-values were calculated by multiplying the Fisher’s exact p-value with the number of tested variants per phenotype. Additional information can be found in S1 File.

Generation of stable cell lines

Variant and wild type miRNA genes were cloned by amplifying the miRNA gene and flanking sequences of at least 130 bp on either side from genomic DNA of patients heterozygous for the variant using sequence specific primers with attachment sites for Gateway cloning. Bead-purified PCR products with homozygous wild type or variant genotype were inserted into the pLenti6/EGFP-V5 destination vector via Gateway cloning as described by Almeida-Souza et al. [48]. Generation and culture of stable SH-SY5Y cells was performed as described previously [49]. Cells were cultured in medium supplemented with 3 μg/ml Blasticidin S (InvivoGen, CA, USA) to select for stable cells. The SH-SY5Y cell line was purchased from Sigma-Aldrich (MO, USA, catalog number: 94030304-1VL). Cell culture media and supplements were purchased from Gibco (Life Technologies, CA, USA). Total RNA was extracted from pelleted cells using the mirVana miRNA Isolation Kit (Ambion, Life Technologies). Additional information can be found in S2 File.

miRNA and EGFP expression analysis

RT-qPCR (reverse transcription–quantitative PCR) was performed using TaqMan microRNA assays or Power SYBR Green PCR Master Mix (Life Technologies) on a ViiA7 Real-Time PCR System (Applied Biosystems, Life Technologies) to assess the expression of mature miRNAs and EGFP in each stable cell line. Data were normalized by geometric averaging of internal control genes [50]. NRQ for miRNAs and EGFP were calculated using the ΔΔCt method based on the approach described by Hellemans et al. [51]. The miRNA NRQ of each clone was normalized for transduction efficiency by division with the EGFP NRQ value. Differential expression was tested using a two-sided unpaired t-test with unequal variance, p < 0.05 was taken as significance threshold. Additional details can be found in S2 File.

Transcriptome analysis

200–400 ng total RNA aliquots of MIR618 cells (clones MIR618WT1-3, MIR618SNP1-3) and MIR204 cells (clones MIR204WT4-5, MIR204SNP4-6) were submitted for analysis with the Human Transcriptome Array 2.0 (Affymetrix, CA, USA) and data normalization using RMA to AROS Applied Biotechnology (Denmark). Prefiltering was performed on the normalized data to exclude transcript clusters lacking a gene symbol or RNA accession number. Differential expression of genes between wild type and mutant cell lines was assessed by SAM v1.0 as implemented in MeV software v4.9.0 [52,53]. A two class unpaired analysis was used (group 1: WT clones, group 2: SNP clones). GSEA was performed using GSEA v2.2.0 software [54,55] using MSigDB v5.0 C5 Gene Ontology molecular function gene sets (c5.mf.v5.0.symbols.gmt) [56]. For all comparisons (group 1: WT clones, group 2: SNP clones), data was collapsed to gene symbols, with max_probe to collapse multiple probes. Gene sets of > 500 and < 15 genes were excluded. 1000 permutations based on gene sets were performed. Gene sets were called significant when FDR < 0.01 and nominal p-value < 0.001.

Small RNA sequencing

200–400 ng total RNA aliquots of MIR618 cells (clones MIR618WT1-3, MIR618SNP1-3) and MIR204 cells (clones MIR204WT4-5, MIR204SNP4-6) were submitted to Biogazelle (Belgium) for small RNA sequencing. Read QC, adaptor trimming, read mapping and annotation to miRBase 20 [22] was performed at Biogazelle for isomiRs and for mature miRNAs (defined as sum of all isomiR reads mapping to a mature miRNA locus). For miRNAs, data was prefiltered with a cutoff of four reads per miRNA per sample. Normalization and differential expression analysis of the miRNA read count data was performed using the DESeq2 Bioconductor package [20]. For DESeq2 analysis of the MIR204 cells the independent filtering option was disabled. miRNAs with Benjamini-Hochberg adjusted p < 0.1 were called significant. For isomiR distribution analysis, isomiR read counts were normalized to the total read count of the mature miRNA (sum of all isomiRs on that locus) in each sample and counts were summed per functional group. IsomiRs were grouped per 5' start sequence: the 5' canonical group contains all isomiRs that have the same 5' end as the canonical miRNA listed in miRBase 20, isomiRs with 5' change are all other isomiRs. The two-tailed unpaired t-test with unequal variance was applied to test for differences in 5’ change isomiR frequencies between genotypes, p < 0.05 was taken as significance threshold.

RNAhybrid prediction

Hybridization of miR-204-5p and miR-204-3p to transcripts was predicted using the RNAhybrid software v2.1.2 [21]. For miRNA targets, up- and downstream flanking regions of 200 nt were used surrounding the hairpin. Prediction was done using an energy cut-off of -20 kcal/mol and a helix constraint from positions 2 to 8 [57].

Data availability

MASTR data has been deposited at the European Genome-phenome Archive, which is hosted by the EBI and the CRG, under accession numberEGAS00001001607. Microarray data is available at Gene Expression Omnibus with accession GSE69675. Small RNA sequencing data are available at Sequence Read Archive with study accession SRP059166.

Supporting Information

S1 File. Supplementary Strategy and Results.

This file contains additional details of the variant screening strategy and additional results (Tables A-E and Figures A-F).

https://doi.org/10.1371/journal.pone.0144428.s001

(DOC)

S2 File. Supplementary Methods.

This file contains additional details of the methods.

https://doi.org/10.1371/journal.pone.0144428.s002

(DOCX)

Acknowledgments

We are particularly grateful to the patients and control individuals that participated in this study. We would like to acknowledge the personnel of the VIB Genetic Service Facility for the genetic analyses. We gratefully acknowledge the assistance of Vicky De Winter in creating lentiviral particles for the generation of stable cells and Rik Hendrickx for IGE sample preparation. We thank Prof. Dr. Vincent Timmerman for the pLenti6/EGFP-V5 destination vector and Dr. Ligia Mateiu for helpful discussions.

Author Contributions

Conceived and designed the experiments: SC MS PDR JDF. Performed the experiments: SC MS. Analyzed the data: SC BS PDR. Contributed reagents/materials/analysis tools: RA AN PDJ SW. Wrote the paper: SC. Interpretation of results: SC MS BS SW PDJ PDR JDF. Diagnosis of individuals: AN RA SW PDJ. Reviewed manuscript: SC MS BS SW AN PDJ RA PDR JDF. Supervised research: JDF.

References

  1. 1. Lee Y, Ahn C, Han J, Choi H, Kim J, Yim J, et al. The nuclear RNase III Drosha initiates microRNA processing. Nature. 2003;425: 415–419. pmid:14508493
  2. 2. Han J, Lee Y, Yeom K-H, Kim Y-K, Jin H, Kim VN. The Drosha-DGCR8 complex in primary microRNA processing. Genes Dev. 2004;18: 3016–3027. pmid:15574589
  3. 3. Yi R, Qin Y, Macara IG, Cullen BR. Exportin-5 mediates the nuclear export of pre-microRNAs and short hairpin RNAs. Genes Dev. 2003;17: 3011–3016. pmid:14681208
  4. 4. Grishok A, Pasquinelli AE, Conte D, Li N, Parrish S, Ha I, et al. Genes and mechanisms related to RNA interference regulate expression of the small temporal RNAs that control C. elegans developmental timing. Cell. 2001;106: 23–34. pmid:11461699
  5. 5. Hutvágner G, McLachlan J, Pasquinelli AE, Bálint E, Tuschl T, Zamore PD. A cellular function for the RNA-interference enzyme Dicer in the maturation of the let-7 small temporal RNA. Science. 2001;293: 834–838. pmid:11452083
  6. 6. Czech B, Hannon GJ. Small RNA sorting: matchmaking for Argonautes. Nat Rev Genet. 2011;12: 19–31. pmid:21116305
  7. 7. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136: 215–233. pmid:19167326
  8. 8. Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? Nat Rev Genet. 2008;9: 102–114. pmid:18197166
  9. 9. Wei Y, Li L, Wang D, Zhang C-Y, Zen K. Importin 8 regulates the transport of mature microRNAs into the cell nucleus. J Biol Chem. 2014;289: 10270–10275. pmid:24596094
  10. 10. Liao J-Y, Ma L-M, Guo Y-H, Zhang Y-C, Zhou H, Shao P, et al. Deep sequencing of human nuclear and cytoplasmic small RNAs reveals an unexpectedly complex subcellular distribution of miRNAs and tRNA 3’ trailers. PLoS ONE. 2010;5: e10563. pmid:20498841
  11. 11. Helwak A, Kudla G, Dudnakova T, Tollervey D. Mapping the human miRNA interactome by CLASH reveals frequent noncanonical binding. Cell. 2013;153: 654–665. pmid:23622248
  12. 12. Pallante P, Visone R, Ferracin M, Ferraro A, Berlingieri MT, Troncone G, et al. MicroRNA deregulation in human thyroid papillary carcinomas. Endocr Relat Cancer. 2006;13: 497–508. pmid:16728577
  13. 13. Feliciano A, Castellvi J, Artero-Castro A, Leal JA, Romagosa C, Hernández-Losa J, et al. miR-125b acts as a tumor suppressor in breast tumorigenesis via its novel direct targets ENPEP, CK2-α, CCNJ, and MEGF9. PLoS ONE. 2013;8: e76247. pmid:24098452
  14. 14. Kan AA, van Erp S, Derijck AAHA, de Wit M, Hessel EVS, O’Duibhir E, et al. Genome-wide microRNA profiling of human temporal lobe epilepsy identifies modulators of the immune response. Cell Mol Life Sci. 2012;69: 3127–3145. pmid:22535415
  15. 15. Feng J, Sun G, Yan J, Noltner K, Li W, Buzin CH, et al. Evidence for X-chromosomal schizophrenia associated with microRNA alterations. PLoS ONE. 2009;4: e6121. pmid:19568434
  16. 16. Mencía A, Modamio-Høybjør S, Redshaw N, Morín M, Mayo-Merino F, Olavarrieta L, et al. Mutations in the seed region of human miR-96 are responsible for nonsyndromic progressive hearing loss. Nat Genet. 2009;41: 609–613. pmid:19363479
  17. 17. Goossens D, Moens LN, Nelis E, Lenaerts A- S, Glassee W, Kalbe A, et al. Simultaneous mutation and copy number variation (CNV) detection by multiplex PCR-based GS-FLX sequencing. Hum Mutat. 2009;30: 472–476. pmid:19058222
  18. 18. Moens LN, De Rijk P, Reumers J, Van den Bossche MJA, Glassee W, De Zutter S, et al. Sequencing of DISC1 pathway genes reveals increased burden of rare missense variants in schizophrenia patients from a northern Swedish population. PLoS ONE. 2011;6: e23450. pmid:21853134
  19. 19. Cammaerts S, Strazisar M, Dierckx J, Del Favero J, De Rijk P. miRVaS: a tool to predict the impact of genetic variants on miRNAs. Nucleic Acids Res. 2015;
  20. 20. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15: 550. pmid:25516281
  21. 21. Rehmsmeier M, Steffen P, Hochsmann M, Giegerich R. Fast and effective prediction of microRNA/target duplexes. RNA. 2004;10: 1507–1517. pmid:15383676
  22. 22. Kozomara A, Griffiths-Jones S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011;39: D152–157. pmid:21037258
  23. 23. Vriens J, Owsianik G, Hofmann T, Philipp SE, Stab J, Chen X, et al. TRPM3 is a nociceptor channel involved in the detection of noxious heat. Neuron. 2011;70: 482–494. pmid:21555074
  24. 24. Chung TKH, Lau TS, Cheung TH, Yim SF, Lo KWK, Siu NSS, et al. Dysregulation of microRNA-204 mediates migration and invasion of endometrial cancer by regulating FOXC1. Int J Cancer. 2012;130: 1036–1045. pmid:21400511
  25. 25. Ying Z, Li Y, Wu J, Zhu X, Yang Y, Tian H, et al. Loss of miR-204 expression enhances glioma migration and stem cell-like phenotype. Cancer Res. 2013;73: 990–999. pmid:23204229
  26. 26. Yin Y, Zhang B, Wang W, Fei B, Quan C, Zhang J, et al. miR-204-5p inhibits proliferation and invasion and enhances chemotherapeutic sensitivity of colorectal cancer cells by downregulating RAB22A. Clin Cancer Res. 2014;20: 6187–6199. pmid:25294901
  27. 27. Wu X, Zeng Y, Wu S, Zhong J, Wang Y, Xu J. MiR-204, down-regulated in retinoblastoma, regulates proliferation and invasion of human retinoblastoma cells by targeting CyclinD2 and MMP-9. FEBS Lett. 2015;589: 645–650. pmid:25647033
  28. 28. Hall DP, Cost NG, Hegde S, Kellner E, Mikhaylova O, Stratton Y, et al. TRPM3 and miR-204 establish a regulatory circuit that controls oncogenic autophagy in clear cell renal cell carcinoma. Cancer Cell. 2014;26: 738–753. pmid:25517751
  29. 29. Conte I, Hadfield KD, Barbato S, Carrella S, Pizzo M, Bhat RS, et al. MiR-204 is responsible for inherited retinal dystrophy associated with ocular coloboma. Proc Natl Acad Sci USA. 2015;
  30. 30. Conte I, Carrella S, Avellino R, Karali M, Marco-Ferreres R, Bovolenta P, et al. miR-204 is required for lens and retinal development via Meis2 targeting. Proc Natl Acad Sci USA. 2010;107: 15491–15496. pmid:20713703
  31. 31. Conte I, Merella S, Garcia-Manteiga JM, Migliore C, Lazarevic D, Carrella S, et al. The combination of transcriptomics and informatics identifies pathways targeted by miR-204 during neurogenesis and axon guidance. Nucleic Acids Res. 2014;42: 7793–7806. pmid:24895435
  32. 32. Paylakhi SH, Moazzeni H, Yazdani S, Rassouli P, Arefian E, Jaberi E, et al. FOXC1 in human trabecular meshwork cells is involved in regulatory pathway that includes miR-204, MEIS2, and ITGβ1. Exp Eye Res. 2013;111: 112–121. pmid:23541832
  33. 33. Xu G, Chen J, Jing G, Shalev A. Thioredoxin-interacting protein regulates insulin transcription through microRNA-204. Nat Med. 2013;19: 1141–1146. pmid:23975026
  34. 34. Hsu S- D, Tseng Y- T, Shrestha S, Lin Y- L, Khaleel A, Chou C- H, et al. miRTarBase update 2014: an information resource for experimentally validated miRNA-target interactions. Nucleic Acids Res. 2014;42: D78–85. pmid:24304892
  35. 35. Huang V, Li L-C. miRNA goes nuclear. RNA Biol. 2012;9: 269–273. pmid:22336708
  36. 36. Matkovich SJ, Hu Y, Dorn GW. Regulation of cardiac microRNAs by cardiac microRNAs. Circ Res. 2013;113: 62–71. pmid:23625950
  37. 37. Tang R, Li L, Zhu D, Hou D, Cao T, Gu H, et al. Mouse miRNA-709 directly regulates miRNA-15a/16-1 biogenesis at the posttranscriptional level in the nucleus: evidence for a microRNA hierarchy system. Cell Res. 2012;22: 504–515. pmid:21862971
  38. 38. Chen P-S, Su J-L, Cha S-T, Tarn W-Y, Wang M-Y, Hsu H-C, et al. miR-107 promotes tumor progression by targeting the let-7 microRNA in mice and humans. J Clin Invest. 2011;121: 3442–3455. pmid:21841313
  39. 39. Butz S, Okamoto M, Südhof TC. A tripartite protein complex with the potential to couple synaptic vesicle exocytosis to cell adhesion in brain. Cell. 1998;94: 773–782. pmid:9753324
  40. 40. Jo K, Derin R, Li M, Bredt DS. Characterization of MALS/Velis-1, -2, and -3: a family of mammalian LIN-7 homologs enriched at brain synapses in association with the postsynaptic density-95/NMDA receptor postsynaptic complex. J Neurosci. 1999;19: 4189–4199. pmid:10341223
  41. 41. Olsen O, Moore KA, Fukata M, Kazuta T, Trinidad JC, Kauer FW, et al. Neurotransmitter release regulated by a MALS-liprin-alpha presynaptic complex. J Cell Biol. 2005;170: 1127–1134. pmid:16186258
  42. 42. Setou M, Nakagawa T, Seog DH, Hirokawa N. Kinesin superfamily motor protein KIF17 and mLin-10 in NMDA receptor-containing vesicle transport. Science. 2000;288: 1796–1802. pmid:10846156
  43. 43. Cheng Q, Zhang X, Xu X, Lu X. MiR-618 inhibits anaplastic thyroid cancer by repressing XIAP in one ATC cell line. Ann Endocrinol (Paris). 2014;75: 187–193.
  44. 44. Fu A, Hoffman AE, Liu R, Jacobs DI, Zheng T, Zhu Y. Targetome profiling and functional genetics implicate miR-618 in lymphomagenesis. Epigenetics. 2014;9: 730–737. pmid:24503492
  45. 45. Nilsson L-G, Bäckman L, Erngrund K, Nyberg L, Adolfsson R, Bucht G, et al. The betula prospective cohort study: Memory, health, and aging. Aging, Neuropsychology, and Cognition. 1997;4: 1–32.
  46. 46. Nilsson L-G, Adolfsson R, Bäckman L, de Frias CM, Molander B, Nyberg L. Betula: A Prospective Cohort Study on Memory, Health and Aging. Aging, Neuropsychology, and Cognition. 2004;11: 134–148.
  47. 47. Reumers J, De Rijk P, Zhao H, Liekens A, Smeets D, Cleary J, et al. Optimized filtering reduces the error rate in detecting genomic variants by short-read sequencing. Nat Biotechnol. 2012;30: 61–68.
  48. 48. Almeida-Souza L, Goethals S, de Winter V, Dierick I, Gallardo R, Van Durme J, et al. Increased monomerization of mutant HSPB1 leads to protein hyperactivity in Charcot-Marie-Tooth neuropathy. J Biol Chem. 2010;285: 12778–12786. pmid:20178975
  49. 49. Strazisar M, Cammaerts S, van der Ven K, Forero DA, Lenaerts A-S, Nordin A, et al. MIR137 variants identified in psychiatric patients affect synaptogenesis and neuronal transmission gene sets. Mol Psychiatry. 2015;20: 472–481. pmid:24888363
  50. 50. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3: RESEARCH0034. pmid:12184808
  51. 51. Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007;8: R19. pmid:17291332
  52. 52. Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A. 2001;98: 5116–5121. pmid:11309499
  53. 53. Saeed AI, Bhagabati NK, Braisted JC, Liang W, Sharov V, Howe EA, et al. TM4 microarray software suite. Methods Enzymol. 2006;411: 134–193. pmid:16939790
  54. 54. Mootha VK, Lindgren CM, Eriksson K- F, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34: 267–273. pmid:12808457
  55. 55. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005;102: 15545–15550. pmid:16199517
  56. 56. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25: 25–29. pmid:10802651
  57. 57. Marín RM, Vanícek J. Efficient use of accessibility in microRNA target prediction. Nucleic Acids Res. 2011;39: 19–29. pmid:20805242