Non-coding RNA dysregulation in the amygdala region of schizophrenia patients contributes to the pathogenesis of the disease

Schizophrenia (SCZ) is a neuropsychiatric disorder with a complex genetic etiology. The redundancy of the gene networks underlying SCZ indicates that many gene combinations have the potential to cause a system dysfunction that can manifest as SCZ or a related neurodevelopmental disorder. Recent studies show that small non-coding microRNA (miRNA) and long non-coding RNA (lncRNA) are important factors in shaping these networks and are dynamically regulated by neuronal activation. We investigated the genome-wide transcription profiles of 46 human amygdala samples obtained from 22 SCZ patients and 24 healthy controls. Using RNA sequencing (RNA-seq), we determined lncRNA expression levels in all samples and generated miRNA profiles for 27 individuals (13 cases and 14 controls). Previous studies have identified differentially expressed miRNAs in SCZ, including miR-132, miR-212, and miR-34a/miR-34c. Here we report differential expression of a novel miRNA, miR1307, in SCZ. Notably, miR1307 maps to a locus previously associated with SCZ through GWAS. Additionally, one lncRNA that was overexpressed in SCZ, AC005009.2, also maps to a region previously associated with SCZ based on GWAS and overlapped SCZ-related genes. The results were replicated in a large independent data set of 254 dorsolateral prefrontal cortex samples from the CommonMind consortium. Taken together, these results suggest that miRNA and lncRNAs are important contributors to the pathogenesis of SCZ.


Introduction
Schizophrenia (SCZ) is a chronic, severe disabling neurodevelopmental disorder that affects people worldwide. While the exact cause of the disorder remains unknown, it is believed to result from an interplay between genetic and environmental factors 1,2 . Recent genetic studies have implicated multiple genomic variants in functional gene networks that individually confer modest risk to the disease rather than high-impact mutations in single genes with large effect sizes [3][4][5][6] . Maintaining optimal network homeostasis requires coordination by an elaborate intracellular network of molecular signal transduction systems, including micro-RNAs (miRNAs), which are involved in the regulation and expression of genetic networks. Therefore, in addition to gene-specific variants that may affect an individual gene or pathway manifesting as SCZ or a related neurodevelopmental disorder 7,8 , miRNAs such as non-coding RNAs that have the potential to dysregulate larger networks of genes have been postulated to be causal factors in the pathogenesis of SCZ 9 . Two groups of miRNAs, including typical non-coding RNAs and long non-coding RNAs (lncRNAs), have been reported to be associated with SCZ as well as other psychiatric diseases 7,10,11 .
The amygdala dysfunction has been implicated in SCZ [12][13][14] ; anatomically it consists of two almond-shaped groups of nuclei located deep and medially within the temporal lobes of the human brain. Functionally, the amygdala is thought to play a primary role in the processing of memory, decision-making, and emotional reactions 15 . Studies in mouse models involving the amygdala suggest that overexpression of miRNAs can protect the animals against anxiety-associated responses 16 and controls their fear response 17 . Research in human subjects also implicates the amygdala region in emotion processing and recognition of facial expressions [18][19][20] . To our knowledge, a genome-wide transcription study of non-coding RNAs in the amygdala of SCZ patients has not been previously reported. In this study, we used nextgeneration sequencing technologies, including RNA sequencing (RNA-seq) and small RNA-seq (SRNA-seq), to assess the expression levels of miRNAs and lncRNAs in postmortem tissue from 46 individuals, including 22 SCZ patients and 24 healthy controls, obtained from the Lieber Institute. Our replication set consisted of 254 dorsolateral prefrontal cortex (DLPFC) tissues obtained from the DLPFC. Data that were generated as part of the Com-monMind Consortium 21 were of European ancestry.

RNA extraction and quality assessment
Total RNA including the miRNA portion from postmortem amygdala tissue was extracted from 100 mg of pulverized tissue with all prep DNA/RNA/miRNA Universal Kit (Qiagen). The RNA portion was purified with RNeasy Mini Spin columns and on-column DNase digestion by RNase-Free DNase Set (Qiagen). Total RNA yield was determined by Qubit RNA BR Assay Kit and Qubit Fluorometer (ThermoFisher Scientific). RNA quality was assessed with high-resolution capillary electrophoresis on an Agilent Bioanalyzer 2100 (Agilent Technologies). Approximately 300 ng of total RNA was applied to an RNA 6000 Nano LabChip without prior heating. An RIN, obtained from entire Agilent electrophoretic trace with the RIN software algorithm, was used for the assessment of RNA quality (scale 1-10, with 1 being the lowest and 10 being the highest RNA quality).

Small RNA library preparation and sequencing
The Illumina TruSeq® Small RNA kit was used to generate the small RNA-seq libraries according to the manufacturers' instructions. Briefly, 3′ (RA3) and 5′ (RA5) adapters were sequentially ligated to Dicer-processed miRNAs. Adaptor-ligated small RNAs were reversetranscribed and the resultant cDNAs PCR-amplified. The PCR products were size-separated on 2% agarose gels and the 22-30 nucleotide bands excised and purified. Purified products were quantified by qPCR before cluster generation on a cBot (Illumina). Twelve samples were multiplexed per lane and sequenced on an Illumina HiSeq2000at 1 × 36 or 1 × 50 bp.
mRNA library preparation and sequencing RNA-seq libraries were constructed using illumina TruSeq RNA sample Prep Kit (RS-122-2001) TruSeq RNA Sample Prep Kit V2, Illumina Inc., San Diego, CA, 92122) following the manufacturers' instructions. The poly-A containing mRNA molecules were purified from 300 to 500 ng DNAse-treated total RNA using oligo (dT) beads. Following the purification, the mRNA was fragmented using divalent cations under elevated temperature (94 degree) for 2 min. Under these conditions, the resultant fragment lengths ranged from 130 to 290 bp with a median length of 185 bp. Reverse transcriptase and random primers were used to generate the first-strand cDNA from the RNA fragments. Second-strand DNA was synthesized using DNA polymerase I and RNaseH. The resulting cDNAs were end-repaired using T4 DNA polymerase, T4 PNK, and Klenow DNA polymerase and a single "A" base added using Klenowexo (3′-5′ exo minus). Illumine PE adapters were ligated using T4 DNA ligase including a barcoding index for multiplexing. These products were then purified and enriched by PCR to create the final cDNA library for high-throughput sequencing on a Highseq2000. The concentration of the RNA-seq libraries was determined by Qubit (Invitrogen, CA) and qPCR, and the quality of the libraries assayed on a Lab ChipGX (Caliper, MA) HT DNA 1 K/12 K/Hisensitivity LabChip. The libraries were multiplexed and loaded on a flow cell for cluster generation on a cBot (Illumina). Illumina Real Time Analysis was used for image analysis and base calling. FASTQ files were generated using the BCL Converter (CASAVA v1.8.2).
RNA-seq/small RNA-seq alignment and differential expression test Following sequencing we obtained an average read depth of 116 M reads per sample, ranging from 34 to 226 M reads across the 46 amygdala samples. Alignments were generated using the Genomic Short-read Nucleotide Alignment Program 22 taking common SNPs (dbSNP137) into account during the alignments. The output SAM files were converted to BAM, sorted by index, and unpaired reads were filtered with SAMtools 23 . Following alignment and filtering an average of 103 M reads were obtained per sample with an average mapping/filtering rate of 89% (Supplementary Table 1). Differential expression tests were performed using the cuffdiff package in cuf-flink2.2.1 24 and the GTF file from GENCODE version 19 25 . The DLPFC RNA-seq data were processed using the same bioinformatics pipeline as previously described for the amygdala data. To reduce algorithm-specific false-positives and to adjust for clinical covariates, we also analyzed the data set using DEseq2 including age, gender, postmortem index (PMI), and manner of death as covariates in the model. Differentially expressed non-coding RNAs were identified from the intersect of the Cuffdiff and DEseq2 26 data sets.
The small RNA-seq data read depths ranged from 4.7 to 72 M with average 20 M reads per sample. Cutadapt1.8.1 was applied to trim the adaptors 27 and the trimmed data were further aligned with Micro Read Fast Alignment Search Tool (MRFast) 28 . The average mapping rate after trimming was 75.1% with an average read depth of 16 M reads per sample.

Clinical sample selection and sample data QC
A detailed description of the clinical sample is provided in Supplementary Table 6. RNA-seq and small RNA data were generated on all 49 individuals listed in Supplementary Table 6. Following primary analysis of the RNAseq data, three individuals (10451, 10678, and 10495) were excluded from further analysis because of outlying expression values. The five most highly differentially expressed genes (CCL3, CCL4, CCL8, CXCL10, and CXCL11) in the analysis were driven by these three controls, with extreme fragments per kilobase million (FPKM) values (~×100 higher than the other samples). The outliers may be related to the cause of death of these individuals, two (10451, 10678) died from "complications of smoke inhalation", and the third (10495) from an accident with multiple head injuries. Of the remaining 46 cases and controls, all were of Caucasian ancestry, 40 were males, and 6 females. Average age at death was 38.4 years. As confounders are an important consideration in postmortem studies, we carried out non-parametric analyses to test for association between age, sex, and PMI with SCZ and gene expression level (FPKM) in the remaining 46 samples. No significant associations were found between FPKM or SCZ and the covariates. Following quality control of the small RNA-seq data, 19 of the 46 individuals with RNA-seq data were excluded from the analysis.
In validation data set DLPFC, RNA was isolated from 50 mg homogenized tissue in Trizol using the RNeasy kit based on the instructional protocol. The mean total RNA yield was 15.3 µg (±5.7). The RNA integrity number (RIN) was determined by fractionating RNA samples on the 6000 Nano chip (Agilent Technologies) on the Agilent 2100 Bioanalyzer. The mean RIN was 7.7 (±0.9), and the mean ratio of 260/280 was 2.0 (±0.02). Processing order was re-randomized prior to ribosomal RNA (rRNA) depletion. Briefly, rRNA was depleted from about 1 µg of total RNA using Ribo-Zero Magnetic Gold kit (Illumina/ Epicenter Cat # MRZG12324) to enrich for polyadenylated coding RNA and non-coding RNA. The sequencing library was prepared using the TruSeq RNA Sample Preparation Kit v2 (RS-122-2001-48 reactions) in batches of 24 samples. A pool of 10 barcoded libraries was layered on a random selection of two of the eight lanes of the Illumina flow cell bridge amplified to~250 million raw clusters. One-hundred base pair paired end reads were obtained on a HiSeq 2500. The sequence data were processed for primary analysis to generate QC values. Samples with a minimum of 50 million mapped reads (~25 million paired end reads) and less than 5% rRNAaligned reads were retained for downstream analysis 21 .

Results
Replication of known SCZ-related miRNA from the small RNA-seq data The small RNA-seq data set included 27 individuals (13 cases and 14 controls); differential expression analysis resulted in the identification of 17 statistically significant differentially expressed miRNAs ( Table 1). Ten of the seventeen are novel discoveries, the remaining seven have previously been implicated in SCZ and/or other psychiatric phenotypes, including two members of miRNA family 34, mir-34c, -34a (p value = 3e−4 and 0.02, respectively) as well as a number of others including mirs-124,-132, -212, -663, and -144.

miRNA based on RNA-seq
While RNA-seq is not designed to systematically detect expression changes of miRNA, some miRNAs that are affiliated with the RNA class are detectable by RNA-seq. Analysis of the GTF file as described in the Materials and methods section resulted in the identification of differential expression of MIR1307 between the SCZ patients and healthy controls. Expression of MIR1307 was approximately threefold higher in controls (p value = 0.0152). The expression difference and the direction of change were confirmed in the small RNA-seq data set, plotting the expression levels in cases and controls confirms that the association is not driven by outliers (Fig.  1a).

LncRNAs based on RNA-seq
Differential analysis of the transcriptome data identified 250 lncRNAs from Gencode version 19 that showed significant expression differences between the amygdala regions of cases and controls (Supplementary Table 2). Among those, two lincRNAs RP11-724N1.1 (p value = 0.014) and RP11-677M14.2 (p value = 0.045) are of particular interest as they map to regions previously associated with SCZ based on genome-wide association study (GWAS) ( Table 2). RP11-724N1.1 resides within the GWAS locus on chr10, previously shown to be associated with SCZ in multiple studies 29,30 . RP11-677M14.2 is located within the SCZ-associated gene, neurogranin  30 , the human homolog of the neuron-specific rat RC3/neurogranin gene, which is a direct target for thyroid hormone in the human brain, and controls its expression.

Non-coding RNAs for SCZ subtypes
We further categorized the SCZ patients based on their clinical subtypes, including disorganized (seven individuals), paranoid (five individuals), and undifferentiated (nine individuals; Supplementary Table 7). The analysis revealed subtype specificity in non-coding RNA expression levels. MIR1307 showed 10-fold higher expression in controls when compared with the undifferentiated subtype in contrast to a two-fold increase when compared with the paranoid subtype. MIR1282 showed significantly higher expression in the paranoid subtype. We also identified subtype specificity in the expression of miRNAs previously reported to be associated with SCZ or other psychiatric diseases. For instance, the miRNA family 34 were only significant in the undifferentiated subtype and mir-144 was only associated  Table 4). For lncRNAs, RP11-677M14.2 was shown to be differentially expressed in the paranoid subtype (p value = 0.0064), whereas RP11-724N1.1 was differentially expressed in undifferentiated SCZ (p value = 0.016; Supplementary Table 5). Finally, some associations were restricted to a particular subtype such as lncRNA AC005009.2, which was not differentially expressed in the combined analysis but is significantly associated with the undifferentiated subtype (p value 0.0016).
Plotting the expression value distributions of the three lncRNAs differentially expressed by subtype (Supplementary Figs. 2a-e) indicates that the results are not driven by outliers. To determine the effects of variables such as age, gender, and manner of death ( Supplementary  Fig. 3), we used expression values of target lncRNAs as dependent variable (y) and age, gender, and mode of death as independent variable (x). For RP11-724N1.1, the p value of linear regression is 0.312, 0.973, and 0.504 for age of death, gender, and manner of death. For RP11-677M14.2 the p values are 0.2735, 0.3565, and 0.0954, and for AC005009.2 the p values are 0.0429, 0.7158, and 0.7331, respectively. Therefore, the expression level differences of subtypes are not driven by the biological parameters at a significance level of p ≤ 0.01. The entire distribution difference p values of expressed lncRNAs based on Wilcoxon rank-sum between undifferentiated, paranoid, and disorganized are 0.056, 0.15, and 0.89. The statistic results indicate that the lncRNA distribution between undifferentiated and paranoid are close to significant difference, and there are certain level differences between paranoid and disorganized, but the differences between undifferentiated and disorganized are insignificant. As a result, it suggests that paranoid SCZ has unique lncRNA patterns compared to other subtypes.
The DESeq2 results of the three lncRNAs were highly concordant with the cuffdiff outputs, expression fold changes were consistent, and almost identical between two tools (Supplementary Table 8). The boxplot of those three lncRNAs showed that the expression differences are not driven by outliers (Figs. 1b-d).

Replication/validation from independent data set
With limited availability of postmortem amygdala tissue from SCZ patients and no previously reported human amygdala RNA-seq studies available in the public domain,  we attempted replication of our findings in DLPFC 21 , an independent data set from the CommonMind Consortium. RNA-seq was carried out in the DLPFC of totally 254 individuals (120 SCZ and 134 controls) of European ancestry. We applied the same bioinformatics analysis pipeline to the DLPFC data as previously described for the amygdala data to validate our lncRNA/miRNA findings. Adequate coverage was obtained for lncRNA AC005009.2 in the replication set but not for the rest three non-coding RNAs. The read depth of DLPFC data set is about half that of the amygdala data at 60~70 M reads vs. 150 M reads. As shown in Table 3, AC005009.2 replicated with a significant p value, and MIR1307 showed a clear trend toward association, with the same direction of effect, in the DLPFC data (Table 3).

Discussion
SCZ is a complex polygenic disease. Multiple risk genes have already been identified through GWAS and sequencing studies 30,31 . To date, no variants, common or rare, conferring large effect sizes have been identified in SCZ. Rather a picture is emerging from the genetic studies implicating the dysregulation of complex gene networks and regulatory mechanisms. Growing evidence indicates that distinct neuronal ncRNAs, particularly miRNAs and lncRNAs, are likely to influence the development of psychiatric diseases, including SCZ. For instance, the lncRNA, Gomafu, is acutely regulated in response to neuronal activation and is involved in SCZ-associated alternative splicing. In this study, we detected altered expression of miRNAs and lncRNAs, well defined in GENCODE annotation, by RNA-seq and SRNA-seq of amygdala tissue samples from 46 subjects including SCZ subjects and controls.
Amygdala dysfunction has been extensively reported in SCZ, yet to our knowledge, gene expression profiling has not been previously reported on human amygdala tissue. Previous studies in mice show that overexpression of miRNA 34 genes in amygdala is associated with anxiety and fear. Overexpression of mir-34c in the amygdala of mice has been shown to induce anxiolytic behavior after challenge 16 ; in our data set, the expression of mir-34c was twofold higher in the SCZ group compared to the control. Mir-34a has been implicated in amygdala-dependent fear memory consolidation in mice through Notch signaling 17 and in our data expression of mir-34a was significantly upregulated in cases. Mir-34a has been shown to be differentially expressed in other brain regions, including the prefrontal cortex 32 , leading to the suggestion that mir-34a levels in blood could serve as a biomarker for SCZ 33,34 . In our study, we identify that the 34 family is only differentially expressed in the individuals with "undifferentiated" subtype of SCZ. Another example is mir-132, a well-known miRNA involved with circadian rhythm and previously associated with multiple psychiatric traits in mouse models 35 and with network regulation in major depressive disorder in human prefrontal cortex tissue 36 . Mir-132 is also differentially expressed in human prefrontal cortex tissue in miRNA-profiling data 32 and is a potential biomarker in human blood 33 . In our study, we show that mir-132 is differentially expressed across all of the clinical SCZ subtypes. Another miRNA of interest, mir-663, is a known regulator of neuronal differentiation 37 . Finally, mir-144 is expressed in the hippocampus of rats and is known to be influenced by both valproate and the mood stabilizer lithium 38 ; mir-144 is found to be only differentially expressed in individuals categorized as "disorganized" SCZ.
In addition to the seven previously reported miRNAs, we further identified 10 miRNAs that will require validation in future human or mouse studies. We also identified over 200 differentially expressed lncRNAs in this  Supplementary Fig. 1a.
CNNM2 encodes a brain-expressed transmembrane protein that is involved in magnesium transport 40 . Loss of function mutations in CNNM2 has recently been reported to be causal of hypomagnesemia, seizures, and mental retardation 41 . RP11-677M14.2 maps to the antisense strand overlapping the NRGN gene ( Supplementary Fig.  1b). The NRGN locus has been shown to be associated with SCZ in the PGC2 meta-analysis (rs55661361, p value = 3.68E−12), as well as in multiple other studies that reported association of functional coding variants in NRGN with SCZ 39, 42-45 . Post-transcriptional regulation of NRGN through overlapping sense and antisense transcripts has previously been reported during cerebral corticogenesis and synapse function in mice 46 . We can confirm through our RNA-seq expression data (Supplementary Table 3) that the expression of NRGN is strikingly higher (FPKM = 1251) in the SCZ group compared to controls (FPKM = 0). Subtype analysis further indicated that RP11-677M14.2 under expression was restricted to cases with a clinical diagnosis of paranoid SCZ. A subtype-specific lncRNA, AC005009.2, has been identified for undifferentiated SCZ individuals; the lncRNA shows twofold higher expressions in control compared to SCZ patients who were diagnosed as the "undifferentiated" subtype. AC005009.2 is also close to GWAS region of PGC2 SCZ meta-analysis (rs12704290, p value = 1.04 × 10 −10 ). Meanwhile, AC005009.2 overlaps the transcript of the metabotropic glutamate receptor 3 gene (GRM3; Supplementary Fig. 5c). GRM3 has been shown to associate with SCZ risks in previous studies [47][48][49] , and AC005009.2 is only differentially expressed in undifferentiated SCZ subtype with p value = 0.0046 and 1.5-fold higher expression in controls.
MIR1307 also maps to a highly significant SCZ GWAS locus (rs11191419, p value = 6.198e−19) in the PGC2 meta-analysis 30 . The associated locus contains multiple genes; however, we only detected differential expression in our RNA-seq data set from MIR1307 (Table 4). The expression data raise the possibility that the GWAS association is being driven by MIR1307. Further work will be required to assess the association between the SCZassociated GWAS variants and MIR1307 expression levels in the brain.
To validate and replicate our findings in the absence of an independent human amygdala RNA-seq data set, we turned to a much larger independent data set of 254 prefrontal cortex tissues. The involvement of the DLPFC in SCZ is well established with multiple reports of association of DLPFC dysfunction correlating with clinical phenotype 50,51 . AC005009.2 replicated with a significant p value and MIR1307 showed a clear trend toward association, with the same direction of effect, in the DLPFC data. The two other lncRNAs, RP11-724N1.1 and RP11-677M14.2, did not reach statistical significance in the DLFPC data, which may either be reflective of tissue specificity or simply due to the reduced coverage in the DLPFC data. The average read coverage in the DLPFC was 70 M reads compared with the 150 M reads from the amygdala data.
Finally, we note that several previously reported SCZassociated miRNAs, including MIR137, 181b, 19, 219, and MIR9, did not replicate in this study. MIR137 and MIR181b were based on GWAS results rather than brain tissues, and MIR181b was reported in Chinese population rather than in Caucasian [52][53][54] . GWAS studies based on blood showed that MIR137 mediated dysregulation as a previously unknown etiologic mechanism in SCZ 52 . Our results showed that MIR137 expression is almost same in SCZ and controls, and MIR181b is 1.2-fold higher than SCZ. MIR19 was reported recently to have abnormal expression in neural progenitor cells from SCZ patientderived induced pluripotent stem cells (iPSCs) 55 ; in our data the expression of MIR19 was 1.2-fold higher in SCZ patients, but it did not reach significance (p value = 0.43). MIR219, a brain-specific miRNA, was reported in both normally developing and SCZ patient iPSC-derived neural stem cells 56 . In our amygdala data, the expression level of MIR219 was 1.4-fold higher in SCZ; however, it did not reach significance with a p value of 0.19. The association of those two miRNA was initially reported in SCZ NPCs but not in brain tissue, which may explain the lack of replication. MIR-9 was reported to be abundantly expressed in control neural progenitor cells (NPCs) and significantly downregulated in a subset of SCZ NPCs 57 . In Cuffdiff pipeline, the mir-9 was marked as region with too many fragments, so no p value was generated. In contrast, DESeq2 reports significant q value of mir-9 (0.02).

Conclusion
Differential expression analysis of small RNAs in SCZ brain identified multiple genes that may be contributing to the pathogenesis of the disease. Intersecting the list of differentially expressed small and lncRNAs with loci identified in previously reported GWAS studies has highlighted two genes, miRNA (mir1307) and lncRNAs (AC005009.2), that are differentially expressed in SCZ and as such may underlie or contribute to the observed GWAS signals at these loci. This study demonstrates the utility of transcriptional profiling of relevant disease tissues in identifying genes underlying GWAS signals.