Identification of the primate-specific gene BTN3A2 as an additional schizophrenia risk gene in the MHC loci.

BACKGROUND: Schizophrenia is a complex mental disorder resulting in poor life quality and high social and economic burden. Despite the fact that genome-wide association studies (GWASs) have successfully identified a number of risk loci for schizophrenia, identifying the causal genes at the risk loci and elucidating their roles in disease pathogenesis remain major challenges.
METHODS: The summary data-based Mendelian randomization analysis (SMR) was used to integrate a large-scale GWAS of schizophrenia with brain expression quantitative trait loci (eQTL) data and brain methylation expression quantitative trait loci (meQTL) data, to identify novel risk gene(s) for schizophrenia. We then analyzed the mRNA expression and methylation statuses of the gene hit BTN3A2 during the early brain development. Electrophysiological analyses of CA1 pyramidal neurons were performed to evaluate the excitatory and inhibitory synaptic activity after overexpression of BTN3A2 in rat hippocampal slices. Cell surface binding assay was used to test the interaction of BTN3A2 and neurexins.
FINDINGS: We identified BTN3A2 as a potential risk gene for schizophrenia. The mRNA expression and methylation data showed that BTN3A2 expression in human brain is highest post-natally. Further electrophysiological analyses of rat hippocampal slices showed that BTN3A2 overexpression specifically suppressed the excitatory synaptic activity onto CA1 pyramidal neurons, most likely through its interaction with the presynaptic adhesion molecule neurexins.
INTERPRETATION: Increased expression of BTN3A2 might confer risk for schizophrenia by altering excitatory synaptic function. Our result constitutes a paradigm for distilling risk gene using an integrative analysis and functional characterization in the post-GWAS era. FUND: This study was supported by the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (XDB02020003 to Y-GY), the National Natural Science Foundation of China (31730037 to Y-GY), and the Bureau of Frontier Sciences and Education, Chinese Academy of Sciences (QYZDJ-SSW-SMC005 to Y-GY).

Evidence before this study

Schizophrenia is a complex, highly heritable and heterogeneous disease with estimated heritability ranging up to 80%. During the past decade, the genome-wide association studies (GWASs) have identified numerous schizophrenia-associated variants and loci in world-wide populations, but most of the risk single-nucleotide polymorphisms (SNPs) or variants identified by GWASs were located in intergenic, intronic or other non-coding regions. To identify the causal genes at the reported risk loci and to elucidate how these risk genes influence the pathogenesis of schizophrenia remains a daunting task.

Added value of this study

By integrating a large-scale schizophrenia GWAS data and brain prefrontal cortex expression quantitative trait loci (eQTL) data, we identified BTN3A2 as a potential schizophrenia risk gene (independent of the C4 gene) in the MHC region. We found that BTN3A2 expression in human brain is highest post-natally. Overexpression of this gene specifically suppressed the excitatory synaptic activity onto CA1 pyramidal neurons and decreased the presynaptic release probability, most likely through its interaction with the presynaptic adhesion molecule neurexins. These results suggested that increased expression of BTN3A2 might confer risk for schizophrenia by altering excitatory synaptic function.

Implications of all the available evidence

We identified BTN3A2 as a risk gene for schizophrenia via integrative analysis of GWAS, eQTL data and brain methylation expression quantitative trait loci (meQTL) data. The role of this primate-specific gene was further characterized by electrophysiological experiments and cell surface binding assays. Our result constituted a paradigm for distilling risk gene using an integrative analysis in the post-GWAS era and provided helpful information for understanding the complex pathogenesis of schizophrenia.

Alt-text: Unlabelled Box

Introduction

Schizophrenia is a complex, highly heritable and heterogeneous disease characterized by episodic psychosis and altered cognitive function [1]. It affects approximately 0.5–1% of the world population and causes significant global burdens. Its morbidity and mortality have been both very high that this disease is widely referred to as the “cancer of mental illness” [2]. So far, the pathophysiology of schizophrenia has not been well elucidated, probably because of the complexity and heterogeneity of the psychopathology and the associated cognitive impairments [1]. An increasing evidence from epidemiology has linked this disorder to prenatal and early postnatal life, and evidence from schizophrenia genetic studies has also suggested that many risk genes of this illness enriched in early neurodevelopmental process such as neuronal differentiation and migration [2]. The ‘neurodevelopmental model’ of schizophrenia posited that perturbations in ‘normal’ brain development, which were associated with molecular changes in the developmental stage and affected by environmental factors, might lead to an altered brain developmental trajectory and consequently cause schizophrenia in early adulthood [2]. However, the underlying neurodevelopmental pathogenesis of schizophrenia has not been well determined.

Due to a high heritability of schizophrenia, many researchers aimed at parsing its pathogenesis through genetic analyses, in particular the genome-wide association studies (GWASs). During the past decade, these GWASs have identified numerous schizophrenia-associated variants and loci in world-wide populations [[3], [4], [5], [6], [7], [8], [9], [10], [11]]. In 2014, the Schizophrenia Working Group of the Psychiatric Genomics Consortium (PGC) carried out a large-scale schizophrenia GWAS (PGC2 GWAS) [10] and reported 108 independent risk loci based on a multi-stage GWAS of >150,000 samples. Through integrating the genetic association signals from the schizophrenia GWAS and expression quantitative trait loci (eQTL), new susceptibility genes have been identified [[12], [13], [14], [15], [16]]. Most recently, Pardiñas et al. [7] carried out another large-scale schizophrenia GWAS (CLOZUK + PGC GWAS) and reported 145 independent risk loci based on a multi-stage GWAS of >100,000 samples. Although GWASs have identified many loci associated with schizophrenia, there has no one-to-one Mendelian mapping between these schizophrenia risk alleles and diagnosis [17]. Instead, schizophrenia is truly complex and appears to result from a myriad of genetic variants each exerting small effects on the overall disease risk, conforming closely to a classical polygenic model [4,18]. Indeed, our recent whole-genome sequencing of monozygotic twins discordant for schizophrenia revealed multiple genetic risk factors for this disease [19]. Approximately 90% of the single-nucleotide polymorphisms (SNPs) or variants identified by GWASs were located in intergenic, intron or other non-coding regions [10]. Elucidating how these risk loci influence the pathogenesis of schizophrenia remains a daunting task. Previous studies showed that most complex disease-associated variants confer risk for the illnesses by acting as eQTL to influence gene expression in cis or in trans [[20], [21], [22], [23]]. Moreover, pathway and integrative analyses showed that the schizophrenia risk genes identified by GWASs were enriched in neurodevelopment-associated pathways like glutamatergic neurotransmission and neuronal calcium signaling pathways, supporting the notation that disturbation of the neurodevelopment processes may play a pivotal role in the pathogenesis of schizophrenia [10,24].

In this study, we aimed to identify genes that are cis-regulated by nearby schizophrenia GWAS risk variants based on a two-stage integrative analysis, followed by functional characterization of the target gene. We identified BTN3A2 as a potential novel schizophrenia risk gene, and found that increased expression of BTN3A2 might confer risk for schizophrenia by altering excitatory synaptic function.

Materials and methods

GWAS associations of schizophrenia

Many GWAS studies of schizophrenia have been conducted in the past decade [3,4,7,10,11]. Because the samples in these GWAS studies were (partially) overlapped [8,11,25] or the dataset were not publically released [[3], [4], [5], [6]], we focused on the result of the recent GWAS study (CLOZUK + PGC GWAS) conducted by Pardiñas et al. [7] in the current integrative study. Briefly, Pardiñas et al. [7] collected genome-wide genotype information for 11,260 schizophrenia cases from the United Kingdom and 24,542 controls from public repositories or through collaborations. They combined these samples with 29,415 cases and 40,101 controls from the PGC2 GWAS [10] and performed meta-analysis, which resulted in a total of 105,318 individuals (including 40,675 schizophrenia cases and 64,643 controls) [7]. More details about sample description, genotyping and statistical analyses could be found in the original paper [7]. We downloaded the summary statistics of this meta-analysis [7] from http://walters.psycm.cf.ac.uk/.

Brain eQTL and meQTL data

The BrainCloud dataset, which was generated using Illumina Human 49 K Oligo array (30,176 probes) and contained 268 human dorsolateral prefrontal cortex (DLPFC) samples collected from fetal (negative ages) through gerontal (80 years old) donors, was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/; GSE30272) [26]. The corresponding SNP data for these subjects (268 individuals), generated using Illumina Infinium II 650 K or Illumina Infinium HD Gemini 1 M BeadChips, were downloaded at the dbGaP (https://www.ncbi.nlm.nih.gov/gap; phs000417.v1.p1) [26]. A total of 654,333 SNPs were successfully genotyped. More detailed information about the BrainCloud dataset can be found in the original paper [26].

For the eQTL analysis, we firstly excluded SNPs with a minor allele frequency (MAF) <0.01 (MAF was estimated on the basis of the 1000 Genomes data [27]) from further analyses. In addition, only probes for which the P-value of the top associated cis-eQTL (<1 Mb away from the probe) was <5 × 10−8 were included in subsequent analyses, as SNPs showing a significant P-value would be strongly associated with the gene expression based on the basic assumption for a Mendelian randomization analysis that the instrumental variable has a strong effect on exposure [28]. In brief, a total of 542,091 SNPs and 4876 probes for eQTL analysis were used. Next, we used Plink (version v1.07) [29] to assess the association between each SNP and each probe. We excluded transcripts with a missing value of >5% of the total samples from the analysis. The following minimum SNP cut-off values were used: per sample call rate should be at least 90%, per SNP call rate should be at least 90%, per SNP MAF should be at least 1%, and lack of significance (P-value >.05) should be observed in the Hardy-Weinberg equilibrium tests. SNP position and gene position were based on the hg19 genome assembly (http://www.ensembl.org/info/data/ftp/index.html).

We used the data of the Psychiatric Encyclopedia of DNA Elements (PsychENCODE) project, which collected prenatal and postnatal post-mortem brain specimens across the developmental continuum to refine and characterize schizophrenia risk loci [30], as an independent brain eQTL dataset for validation. We downloaded the eQTL results with a false discovery rate < 0.05 and a filter requiring genes to have an expression >0.1 FPKM (fragments per kilobase per million mapped fragments) in at least 10 samples from the PsychENCODE Integrative Analysis website (http://resource.psychencode.org/) [31]. The PsychENCODE Integrative Analysis [31] contained adult brain prefrontal cortex data of 1387 individuals from the PsychENCODE [30] and the Genotype-Tissue Expression (GTEx, https://www.gtexportal.org) data [32]. Detailed information about the data collection and analysis process can be found in the original study [31].

We also retrieved the release v.7 summary statistics from the eQTL data of 14 brain-related tissues of the GTEx project [32]. The GTEx project is a resource database and associated tissue bank for the scientific community to study the relationship between genetic variation and gene expression in human tissues [32]. Currently, the release of GTEx (release v.7) contains a total of 11,688 samples from 53 tissues of 714 donors within an age range from 20 years to 79 years [32]. In brief, the RNA expression data was collected using Illumina TrueSeq RNA sequencing (non-stranded, ployA+ selection) and Affymetrix Human Gene 1.1 ST Expression Array. Whole genome/exome sequencing and Illumina SNP Array were then used to get the genotype data. More detailed summary and data processing method of the dataset can be found in the GTEx website (https://www.gtexportal.org) [32].

The binary format methylation quantitative trait loci (meQTLs) datasets were downloaded from the summary data-based Mendelian randomization analysis (SMR) website (http://cnsgenomics.com/software/smr) [28]. Briefly, Qi et al. [33] conducted a meta-analysis of meQTL data from ROSMAP [34], Hannon et al. [35] and Jaffe et al. [36]. These meQTL data in binary format were stored at the SMR website [28]. The detailed information about the sample size and method can be found in the original publication of meQTL meta-analysis [33] and the SMR website [28].

SMR integrative analysis

The basic assumption of SMR [28] is if the expression level of a gene is influenced by an eQTL, then there will be differences in gene expression levels among individuals carrying different genotypes of the genetic variant. Then, if the expression level of the gene has an effect on a trait, we will observe differences in phenotype among the different genotype groups: that is, the genetic variant will also show an effect on the trait [28]. The SMR employs the Mendelian randomization to test for pleiotropic associations between gene expression and complex traits using eQTL and trait GWAS summary data [28]. This method can be used to integrate data from multiple-omics studies, such as GWAS data and meQTL data [28].

We integrated the eQTL data [26,32,33] and the meQTL data [33] with the GWAS summary data [7] to perform the SMR integrative analysis [28]. To prepare these data, a sparse BESD file was made by SMR (version 0.631) [28] using the –make-besd argument, which provides an efficient way to store the eQTL summary data in binary format. The sparse BESD file stored the data for eQTLs within 2 Mb of a probe in either direction, as well as SNPs within 1 Mb of any eQTLs with P-value <1 × 10−5 in either direction in the rest of the genome. For GWAS data [7], we removed the SNP with a MAF <0.01 (MAF was estimated based on the 1000 Genomes data [27]). SMR (version 0.631) [28] was then used to integrate the eQTL (or meQTL) data (BESD format) and GWAS data [7] using the following settings: –bfile (individual-level SNP genotype data in PLINK binary format), −-gwas-summary (GWAS summary data) and –beqtl-sumary (BESD format eQTL or meQTL data).

Sherlock integrative analysis

Sherlock (http://sherlock.ucsf.edu/submit.html) is a powerful integrative method and computes the Bayes factor for each gene [37]. The basic assumption of Sherlock is that the expression level of a specific gene(s) may influence the risk of a disease. For a given gene, there may be many variants in the genome affecting its expression (expression SNPs, eSNPs). A change of genotype at any of these eSNPs would lead to mRNA expression change of the gene, which could in turn affect the disease risk. Sherlock computes the LBF score for each eSNP of the gene being tested to represent how strongly the SNP supports a functional role of the gene [37]. The total LBF score of a gene evaluates evidence supporting that the gene is associated versus not associated with the disease. The detailed modeling assumption and calculation theory are described in the original study [37].

We firstly used the standalone version of Sherlock [37] to integrate the BrainCloud eQTL data [26] with GWAS data [7] under the following settings: N_expr = 268 (sample size of the BrainCloud eQTL data), is_pheno_binary = 1 (whether the phenotypic trait is binary (1) or quantitative (0)), N_pheno = 105,318 (40,675 + 64,643 = 105,318), K = 0.01 (the disease prevalence). Other settings were set as the default arguments. The integration of GWAS data [7] with the GTEx eQTL data [32] were conducted using Sherlock webserver (http://sherlock.ucsf.edu/submit.html). We uploaded the GWAS data [7] and selected GTEx eQTL (release v.7) [32] data related to brain tissues as the eQTL data sets.

Gene expression and methylation data

We explored the gene expression pattern from Human Brain Transcriptome (HBT, http://hbatlas.org/) [38] and the PsychENCODE Human Brain Development website (http://development.psychencode.org) [39]. HBT is a public database containing transcriptome data and associated metadata obtained from human brain tissues at various development stages [38]. A total of 1340 tissues collected from 57 postmortem human brains from 16 brain regions were sampled, and the exon-level transcriptome data was generated using Affymetrix GeneChip Human Exon 1.0 ST Arrays. The detailed information about the sample origin, sample size and data processing can be found in the original publication [38]. Recently, the PsychENCODE generated transcriptomic profiling data (mRNA-seq) of 607 histologically verified, high-quality tissue samples from 16 anatomical brain regions that were dissected from 41 brains with age range from 8 postconceptional weeks to 40 postnatal years [39]. These data were systematically analyzed and made queryable in the PsychENCODE website (http://development.psychencode.org) [39]. The detailed information about sample origin, sample size and data processing can be found in the original publication [39]. Jaffe et al. [36] generated DNA methylation data of the DLPFC brain tissues from 335 non-psychiatric controls across the lifespan and 191 patients with schizophrenia by using the Illumina HumanMethylation450 microarray, and they identified widespread changes in the transition from prenatal to postnatal life [36]. The sample size was expanded to 450 psychiatric disorder-free samples and 225 schizophrenia samples based on the downloaded data (GEO accession number GSE74193). Detailed information about the sample source, sample size and data processing can be found in the original publication [36] and the GEO dataset GSE74193.

CommonMind Consortium data and human induced pluripotent stem cells (hiPSC) expression data

The normalized DLPFC gene expression profile of schizophrenia and control subjects were downloaded from Synapse (https://www.synapse.org, syn5607581) from the CommonMind Consortium (CMC) project [40] with permissions. This dataset contained mRNA expression data of the post-mortem human brain specimens from three brain banks (the Icahn School of Medicine at Mount Sinai, the University of Pennsylvania and the University of Pittsburg). Briefly, the total RNA was isolated from 258 schizophrenia cases and 279 control subjects, and was sequenced by Illumina HiSeq 2500 after quality control by RNA integrity number. The cleaned sequencing reads were mapped to human reference genome hg19 (http://www.ensembl.org/info/data/ftp/index.html) and gene expression levels were quantified using log(CPM) (read counts per million total reads) [40].

hiPSC gene expression profiles [41] were downloaded at GEO database (GSE25673). This dataset contains the neurons differentiated from hiPSC that were derived from the primary human fibroblasts from 4 schizophrenia patients and 4 control subjects. Most of these hiPSC neurons were presumably glutamatergic and the others were GABAergic and dopaminergic [41]. Gene expression analyses were performed on 6-week-old hiPSC neurons using Affymetrix Human Gene 1.0ST arrays. Three independent neural differentiations from fibroblasts for each donor were compared. More details about the expression of hiPSC neurons can be found in the original paper [41].

Antibodies and cell culture

The following antibodies and cells were used in this study: mouse monoclonal anti-Flag (Eno-Gene, E12-001), peroxidase-conjugated anti-mouse antibody (KPL, 474–1086), Alexa Fluor 488 conjugated donkey anti-mouse secondary antibody (Invitrogen, A21202). HEK293T cells were obtained from the Kunming Cell Bank of Kunming Institute of Zoology. Cells were maintained in DMEM medium (Gibco-BRL, 11965–092) supplemented with 10% fetal bovine serum (Gibco-BRL, 10099–141) at 37 °C in a humidified atmosphere incubator with 5% CO2.

Plasmid construction and transfection

The full length human BTN3A2 (BTN3A2-L, GenBank accession number NM_001197247.2) and a truncated BTN3A2 lacking transmembrane domain (BTN3A2-S, XM_005248831.4) were cloned into Flag-tagged pCMV-3Tag-8 with HindIII and XhoI for constructing overexpression vectors pCMV-BTN3A2-L-Flag and pCMV-BTN3A2-S-Flag, respectively. Mouse Nlgn1 (NM_138666.4), Nlgn3 (NM_172932.4), Nrxn1β (NM_001346959.1), Nrxn1βA4(−) (NM_001346960.1) and Nrxn3β (NM_001252074.2) were also cloned into pCMV-DsRed2, respectively. The cDNA of full length BTN3A2 was also sub-cloned into pCAGGS-IRES-GFP vector for biolistic transfection. All the primers used for constructing plasmids were listed in Table S1. Each construct was verified by direct sequencing.

HEK239T cells were transfected using Lipofectamine™ 3000 (Invitrogen L3000008) following the manufacture's protocol. In brief, cells were cultured in 6-well plates and allowed to grow to 70% confluence. Then cells in each well were transfected with a total volume of around 270 μL mixture containing 2.5 μg DNA, 250 μL Opti-MEM medium (Gibco-BRL, 31985–070) and 12.5 μL Lipofectamine™ 3000. The cells were incubated with the transfection mixture for 6 h, and then the medium was changed into growth medium for 42 h before the harvest.

Electrophysiology in rat hippocampal slice cultures

Electrophysiology in rat slice cultures was conducted as described in our previous studies [42,43]. Briefly, organotypic hippocampal slice cultures were made from postnatal 6 days – postnatal 8 days rats. The vectors were transfected on days in vitro (DIV) 2 during culturing by a Helios Gene Gun (Bio-Rad, America) with 1 μm DNA-coated gold particles. Slices were maintained at 34 °C with media changes every other day. On DIV 8, the responses of pyramidal neurons, which were identified by morphology and location, in area CA1 were simultaneously recorded from a fluorescent transfected neuron and a neighboring untransfected control one. Dual whole-cell recordings measuring evoked excitatory or inhibitory postsynaptic responses used artificial cerebrospinal fluid (ACSF) bubbled with 95% O2/5% CO2 consisting of 119 μM NaCl, 2.5 μM KCl, 4 μM CaCl2, 4 μM MgSO4, 1 μM NaH2PO4, 26.2 μM NaHCO3, and 11 μM glucose. A total of 100 μM picrotoxin was added to block inhibitory currents when measuring AMPA and NMDA excitatory postsynaptic responses (EPSCs) or 10 μM NBQX/50 μM APV when measuring inhibitory postsynaptic responses. In addition, 2-chloroadenosine (4 μM) was used to control epileptiform activity. A bipolar stimulation electrode was placed in stratum radiatum, and responses were evoked at 0.2 Hz. Peak α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor (AMPAR) currents were recorded at −70 mV, and N-methyl-d-aspartate receptor (NMDAR) current amplitudes 100 ms following the stimuli were recorded at +40 mV. Paired-pulse ratio was determined by delivering two stimuli 40 ms apart and dividing the peak response to stimulus 2 by the peak response to stimulus 1. Peak GABA receptor-mediated currents were recorded at 0 mV. All these data were analyzed by custom software (IGOR Pro).

Cell surface binding assay

Cell supernatant containing the truncated isoform of BTN3A2 tagged by Flag (BTN3A-S-Flag) was produced from the HEK293T cells after transfection with expression vector pCMV-BTN3A2-S-Flag (Table S1) for 48 h. The supernatant was added to the transfected HEK293T cells overexpressing DsRed2, Nlgn1-DsRed2, Nlgn3-DsRed2, Nrxn1β-DsRed2, or Nrxn3β-DsRed2 (Table S1), respectively. After an incubation for 4 h, the medium was removed and cells were fixed with 4% paraformaldehyde for 15 min. Cells were then washed for three times with phosphate-buffered saline (PBS) and incubated at room temperature for 30 min in a blocking solution containing 5% bovine serum albumin in PBS. Mouse anti-Flag antibody was then added (Eno-Gene, E12–001, 1:500 ratio) and incubated for 2 h. The slides were washed three times with PBS and incubated with anti-mouse Alexa Fluor 488 fluorescent secondary antibody (Invitrogen, A21202, 1:500 ratio) for 1 h at room temperature to label BTN3A2. The slides were visualized under an Olympus FluoViewTM 1000 confocal microscope (Olympus, America).

Statistical analysis

The mRNA expression levels of different genes in CMC DLPFC [17] and hiPSC neurons [41] were compared using the Student's t-test using R software (https://www.r-project.org/). For the electrophysiology experiment, all the statistical analysis was compared with the respective control neurons using two-tailed Wilcoxon signed-rank sum test. A P-value <.05 was considered to be statistically significant.

Data availability

The datasets used and/or re-analyzed in this study are available at the original sources [7,17,26,[31], [32], [33],38,39,41] and could be found in our webserver SZDB (www.szdb.org/integrative1.php) [12].

Ethics statement

All animal experiments were performed in accordance with established protocols approved by the Institutional Review Board of Kunming Institute of Zoology, Chinese Academy of Sciences.

Results

Identification of BTN3A2 as a schizophrenia risk gene by SMR integrative analysis

The overall study design and relevant rationale were shown in Fig. 1. The genotype and brain DLPFC expression data of 268 individuals were retrieved from BrainCloud [26], which contained 542,091 SNPs (each SNP with a MAF ≥ 0.01) and 4876 probes (each probe with a PeQTL < 5 × 10−8), and were used for the following analyses. As shown in Table 1, BTN3A2, a primate-specific gene [44,45], which was located in the extended major histocompatibility complex (MHC) region, had the most significant association with schizophrenia (PSMR = 2.59 × 10−11, Bonferroni corrected P-value <.05). SNP rs1979, which was located in the 3’UTR region of BTN3A2, significantly affected the mRNA expression of BTN3A2 (PeQLT = 1.44 × 10−17, Fig. S1) [26] and was significantly associated with schizophrenia (P-value = 1.15 × 10−26) [7]. Another previously uncharacterized non-coding RNA (ncRNA) LOC100131289 also survived the SMR association test (PSMR = 4.11 × 10−6, Bonferroni corrected P-value <.05; Table 1). So far, there are no reports for potential function of this ncRNA according to the PubMed search with keyword “LOC100131289”. This ncRNA has a relatively higher expression in digestive system, spleen, ovary / testis and skin, than that of brain (https://www.ncbi.nlm.nih.gov/gene/?term=LOC100131289). Further study should be carried out to discern its potential role in schizophrenia.

Fig. 1

Workflow of integrative analyses and functional assay in this study. The SMR [28] was used to integrate BrainCloud eQTL data [26] and CLOZUK + PGC GWAS data [7]. The Sherlock [37] was used to integrate BrainCloud [26] and GTEx brain tissue eQTL data [32] with CLOZUK + PGC GWAS data [7] to validate the result of SMR. Brain meQTL data [36] was also used to identify schizophrenia risk genes. Next, brain expression data [26,38,39,76] and methylation data [36] were explored to check the expression and methylation pattern during the prenatal and postnatal neurodevelopment. Data from GEO (GSE25673) [41] and CMC (syn5607581) [17] were used to plot the BTN3A2 mRNA expression pattern between normal controls and schizophrenia cases. The electrophysiology assays in rat brain slice and cell surface binding assays were conducted to characterize the potential function of this gene.

Table 1

Integrative analysis (SMR) of schizophrenia GWAS (CLOZUK + PGC GWAS) and BrainCloud eQTL reveals BTN3A2 as a schizophrenia risk gene.

Probe	Gene	Supporting SNP (trans or cis)	PGWASa	PeQTLb	PSMRc	PHEIDId	Padje
hHR012536	BTN3A2	rs1979 (cis)	1.15 × 10−26	1.44 × 10−17	2.59 × 10−11	2.64 × 10−4	< 0.05
hHR010906	LOC100131289	rs175597 (cis)	1.19 × 10−30	5.00 × 10−7	4.11 × 10−6	1.54 × 10−2	< 0.05
hHC024109	GPR75-ASB3	rs3770407 (cis)	6.10 × 10−5	4.53 × 10−25	1.84 × 10−4	6.50 × 10−2	NS
hHC030929	THAP5	rs40919 (cis)	1.02 × 10−4	2.65 × 10−37	1.98 × 10−4	3.73 × 10−3	NS
hHC007598	NDUFAF7	rs1056021 (cis)	1.59 × 10−4	5.58 × 10−46	2.58 × 10−4	3.72 × 10−2	NS
hHC028708	THAP5	rs40926 (cis)	9.95 × 10−5	2.14 × 10−23	2.89 × 10−4	1.98 × 10−3	NS
hHR014285	LINC00339	rs2143103 (cis)	7.04 × 10−6	3.69 × 10−9	3.51 × 10−4	NA	NS
hHR021334	BOP1	rs3802232 (cis)	8.51 × 10−7	2.74 × 10−7	3.77 × 10−4	7.05 × 10−3	NS
hHA037802	SRA1	rs10042299 (cis)	8.03 × 10−5	1.63 × 10−15	4.04 × 10−4	2.26 × 10−2	NS

P-values for eQTL SNPs from the CLOZUK + PGC GWAS data [7].

P-values for eQTL SNPs from the BrainCloud dataset [26].

P-values from summary data-based Mendelian randomization (SMR) analysis.

P-values from heterogeneity in dependent instruments (HEIDI) test. HEIDI uses multiple SNPs in a cis-eQTL region to distinguish pleiotropy from linkage. NA, not available.

We used the Bonferroni correction to account for multiple testing, which resulted in a genome-wide significance level of P = 1.03 × 10−5 (= 0.05/4876, 4876 probes were retained for SMR analysis). NS, not significant.

Validation of BTN3A2 as a schizophrenia risk gene

We conducted a different integrative analysis (Sherlcok [37]) so as to use GWAS data (CLOZUK + PGC GWAS) [7] and eQTL data (BrainCloud eQTL) [26] to validate the results of the SMR integrative analysis. The Sherlock method [37], which computes the Bayes factor for each gene, reasoned that for a disease-associated gene, any genetic variation that perturbs its expression is also likely to influence the disease risk. We could well validate the result that BTN3A2 was a schizophrenia risk gene. By integrating CLOZUK + PGC GWAS data [7] with BrainCloud eQTL data [26] using Sherlock [37], we confirmed BTN3A2 as a schizophrenia risk gene with P-value = 1.68 × 10−7 and the logarithm of Bayes factors (LBF) value of 9.38 (Table S2), which means that the posterior probability of the association between BTN3A2 and schizophrenia was exp. (9.38) = 11,849 times more likely than the opposite hypothesis (no association). By integrating CLOZUK + PGC GWAS data [7] with eQTL data of the PsychENCODE brain PFC [31] and the GTEx 14 brain tissues [32], respectively, we also found that BTN3A2 was significantly associated with schizophrenia (Tables S3 and S4).

The Sherlock [37] analysis of CLOZUK + PGC GWAS data [7] and GTEx eQTL data [32] also confirmed that BTN3A2 was a schizophrenia risk gene (all P-value <1.00 × 10−5), with all LBF larger than 7 in all studied brain regions (Tables S5). All the integrative analysis data could be found in SZDB (www.szdb.org/integrative1.php) [12]. In addition, we investigated the eQTL of rs1979 in Braineac [46] and found that the risk allele G of rs1979 could significantly elevate the mRNA expression of BTN3A2 (Fig. S3), which was consistent with our observations in BrainCloud [26] (Fig. S1) and GTEx [32] (Fig. S2) datasets.

Taken together, we validated that BTN3A2 was a schizophrenia risk gene using independent integrative approaches (SMR [28] and Sherlock [37]), GWAS (CLOZUK + PGC GWAS) [7] and brain eQTLs (BrainCloud [26], PsychENCODE [31] and GTEx [32]) (Fig. 1). More importantly, we consistently detected the elevated mRNA expression levels of BTN3A2 in brain tissues from carriers of the risk allele G of rs1979, suggesting that increased expression of BTN3A2 was a potential risk factor for schizophrenia.

BTN3A2 was independent of the C4 signal in the MHC region

Of note, integrating GWAS data [7] with GTEx eQTL data [32] also identified another MHC region gene C4A as a schizophrenia risk gene (Tables S6), which has been previously reported to mediate synapse elimination during postnatal development and therefore affects the risk of schizophrenia [47]. The MHC region is the most complex region in human genome because of its unintelligible linkage disequilibrium (LD), and this region exhibited the strongest genetic association with schizophrenia [7,10]. We therefore tested whether BTN3A2 was an independent signal or was just a hitchhiking effect to the C4A gene by analyzing the LD relationship between the respective SNPs supporting the involvement of BTN3A2 and C4A in the risk of schizophrenia. Interestingly, these SNPs were distributed in two evident blocks, with BTN3A2 supporting SNPs in one block and C4A supporting SNPs in the other block (Fig. S4). This result suggested that BTN3A2 might be an independent schizophrenia risk gene in the MHC region.

Integrative analysis of methylation and GWAS data indicated BTN3A2 as a schizophrenia risk gene

Through integrative analysis, we found multiple GWAS SNPs conferring the risk of schizophrenia by affecting the mRNA expression of BTN3A2 (Table S4). As elaborated previously, noncoding risk variants might leverage the expression of nearby genes by affecting the binding of trans-acting factors or changing the epigenetic conditions of a cis region, like histone and/or chromatin modification, or DNA methylation [2]. Based on these hypotheses and observations, we used SMR [28] to integrate schizophrenia CLOZUK + PGC GWAS data [7] with brain-mMeta meQTL data [33], and observed significant association between BTN3A2 and schizophrenia (Table 2). Taken together, our integrative analyses combining GWAS data [7] with gene expression [26,31,32] and methylation profiles [33] indicated that GWAS SNPs likely affected the mRNA expression levels of BTN3A2 through modulating the methylation statuses of its adjacent sites.

Table 2

Integrative analysis (SMR) of schizophrenia GWAS (CLOZUK + PGC GWAS) and ROSMAP Hannon Jaffe brain meta meQTL reveals BTN3A2 as a schizophrenia risk gene.

Gene	Probe	Supporting SNP (cis or trans)	PGWASa	PeQTLb	PSMR	PHEIDI	Padj
BTN3A2	cg23465465	rs1977 (cis)	6.63 × 10−27	1.76 × 10−21	1.06 × 10−12	1.16 × 10−2	<0.05
BTN3A2	cg14345882	rs9393714 (cis)	6.52 × 10−27	<1 × 10−20	4.49 × 10−25	3.06 × 10−5	<0.05
BTN3A2	cg02045355	rs3734536 (cis)	9.20 × 10−15	<1 × 10−20	3.03 × 10−14	1.04 × 10−6	<0.05
BTN3A2	cg10795676	rs34565965 (cis)	2.21 × 10−13	2.83 × 10−10	1.73 × 10−6	2.70 × 10−6	NS

P-values for eQTL SNPs from the CLOZUK + PGC GWAS data [7].

P-values for eQTL SNPs from the ROSMAP Hannon Jaffe brain meta meQTL methylation data [33].

Inverse pattern of methylation level and BTN3A2 mRNA expression level during fetal brain development

Transcriptomic and epigenomic profilings associated with schizophrenia might affect the early neurodevelopmental trajectory [39]. A recent report from the PsychENCODE Consortium revealed that schizophrenia risk genes were enriched in a co-expression module which was significantly associated with neurodevelopment, and expression pattern was significantly changed during the prenatal and postnatal neurodevelopment [39]. Therefore, we explored the dynamic pattern of mRNA expression and methylation conditions during the fetal neurodevelopment. We found that BTN3A2 was upregulated during the fetal brain development (Fig. 2A), but the methylation levels of the above described schizophrenia-associated loci (Table 2) showed an opposite trend (except for cg10795676, Fig. 2B). These results further suggested that GWAS risk SNPs might affect the expression of BTN3A2 through altering the methylation profiles.

Fig. 2

BTN3A2 expression in human brain is highest post-natally. (A) BTN3A2 expression pattern in human brain according to PsychENCODE [39], Human Brain Transcriptome [38], BrianCloud [26] and Brainspan [76] datasets. (B) Level of methylation in four sites identified by meQTL integrative analysis based on the data reported by Jaffe et al. [36].

BTN3A2 was differentially expressed between schizophrenia cases and controls

Post-mortem gene-expression studies of schizophrenia patients and controls suggested subtle abnormalities in multiple brain regions (including the prefrontal and temporal cortex, hippocampus, etc.) and several specific cell types [17,48,49]. Given the fact that integrative analyses identified BTN3A2 as a risk gene for schizophrenia and the risk allele G of rs1979 was significantly associated with a higher mRNA expression of BTN3A2, we speculated that this gene might be up-regulated in schizophrenia cases. To test this hypothesis, we downloaded the gene expression data from CMC [17] which sequenced RNA from DLPFC of both schizophrenia patients (N = 258) and control subjects (N = 279) [17]. We found that BTN3A2 was significantly up-regulated (P-value = .018) in schizophrenia cases compared with control subjects (Fig. 3A).

Fig. 3

BTN3A2 is differentially expressed between schizophrenia cases and controls. The mRNA level of BTN3A2 is significantly up-regulated in (A) post-mortem dorsolateral prefrontal cortex tissues and (B) iPSC neurons from schizophrenia patients compared to that of healthy controls. Data were taken from CMC [17] and GSE25673 [41], respectively. The P-values were calculated by using unpaired Student's t-test.

To further validate this result, we downloaded the gene expression data obtained from neurons differentiated from the human induced pluripotent stem cells [41]. This dataset contained microarray data of neurons that were differentiated from reprogrammed fibroblasts from schizophrenia patients and control individuals. The mRNA expression levels of BTN3A2 in neurons were compared between cells derived from schizophrenia cases and control subjects, and we again observed significant increase of the mRNA levels of this gene in schizophrenia cases (Fig. 3B). Overall, these results were consistent with our findings of the integrative analysis, and suggested that a higher expression level of BTN3A2 might increase the risk of schizophrenia.

BTN3A2 specifically regulated excitatory synaptic transmission

Our analyses strongly suggested that BTN3A2 plays essential roles in schizophrenia pathogenesis. We therefore further explored the function of BTN3A2 to provide insights into its role in this illness. BTN3A2 is an immunoglobulin (Ig) superfamily receptor protein of the butyrophilin and butyrophilin-like families, and is known to play critical immunomodulatory roles [45]. Moreover, Ig superfamily adhesion molecules are among the most abundant trans-synaptic signaling molecules and exert diverse functions during synaptic physiological and pathological processes [[50], [51], [52]]. To study the function of BTN3A2 during synaptic transmission, we used cultured rat hippocampal slice as the study system and expressed BTN3A2 exogenously in CA1 neurons through biolistic transfection. Then dual whole-cell recordings were applied to measure the evoked excitatory or inhibitory postsynaptic responses. We found that overexpression of BTN3A2 decreased both AMPAR and NMDAR-mediated excitatory synaptic transmission compared with the control neurons (Fig. 4A and B), while no effect was seen on the ratio of AMPA and NMDA EPSCs (Fig. 4C). In addition, the paired-pulse ratio (Fig. 4D), a parameter for presynaptic release probability, was increased in BTN3A2-overexpressed neurons. These results suggested that BTN3A2 regulated excitatory synaptic transmission through decreasing the presynaptic glutamate release. In contrast, BTN3A2 had no effect on GABA receptor-mediated inhibitory postsynaptic transmission (Fig. 4E). All these results suggested that BTN3A2 was specifically involved in excitatory synaptic maintenance and its dysregulated expression might lead to malfunction of excitatory/inhibitory balance.

Fig. 4

BTN3A2-mediated regulation is specific to excitatory synapses. Overexpression of BTN3A2 decreases both AMPAR (A) and NMDAR (B)-mediated synaptic transmissions compared with the control neurons. Overexpression of BTN3A2 increases the paired-pulse ratio, a parameter for presynaptic release probability (C), but does not change the ratio of AMPAR and NMADR-mediated EPSCs related to neighboring wild-type neurons (D), and has no effect on GABA receptor-mediated inhibitory postsynaptic transmission (E). Differences between neurons overexpressing BTN3A2 and respective control neurons are compared using two-tailed Wilcoxon signed-rank sum test. *P-value <.05, **P-value <.01, ***P-value <.001. Bars represent mean ± SEM.

BTN3A2 interacted with pre-synaptic neurexins

Neurexins (Nrxns) are a family of presynaptic transmembrane proteins that function as the organizer of synapses [53]. They play important roles in presynaptic release regulation and have been implicated in schizophrenia by human genetic studies [[54], [55], [56], [57]]. Therefore, we hypothesized that BTN3A2 might interact with presynaptic neurexins and thereby affect the excitatory synaptic maintenance and excitatory/inhibitory balance [53]. To test this hypothesis, we used cell surface binding assay to test whether BTN3A2 could interact with neurexins. The Flag-tagged and truncated BTN3A2 (BTN3A2-S-Flag), which lacked the transmembrane domain (Fig. 5A), was successfully expressed in HEK293T cells and its secretion to extracellular supernatant was confirmed by Western blot (Fig. 5B). We then overexpressed DsRed2-tagged Nrxn1β or Nrxn3β in HEK293T cells for 48 h, and replaced the extracellular medium with culture supernatant containing secreted BTN3A2-S-Flag for another 4 h. Immunostaining was performed to examine the surface expression of Nrxns and BTN3A2. We found that BTN3A2 was co-localized with either Nrxn1β or Nrxn3β at the cell surface (Fig. 5C), suggesting that they were directly interacted through the Nrxns extracellular domain. However, the alternative splicing at the 4th canonical sites of Nrxn1β (Nrxn1βΑ4 (−)) that is critical for its binding to its postsynaptic ligand, was not required for this interaction as Nrxn1βΑ4 (−) still significantly bound with BTN3A2 (Fig. 5C). The same strategy was used to determine the interaction between BTN3A2 and postsynaptic neuroligins (Nlgns). We found that neither Nlgn1 nor Nlgn3 interacted with BTN3A2 (Fig. 5C). These results suggest that BTN3A2 specifically interacts with presynaptic Nrxns. Based on these observations, we propose a model that BTN3A2 functions as a postsynaptic ligand for neurexins to regulate presynaptic glutamate release, thereby maintaining the excitatory/inhibitory balance of our brain. Abnormality of BTN3A2 expression might cause an imbalance of synaptic transmission and lead to schizophrenia (Fig. 5D).

Fig. 5

Interaction of BTN3A2 with presynaptic neurexins and a working model of the role of BTN3A2 in synapses. (A) Schematic profile of the two isoforms of BTN3A2. The short isoform BTN3A2-S has a deletion of the transmembrane domain of the full length BTN3A2 (BTN3A2-L). (B) The overexpressed BTN3A2-S in HEK293T cells can be secreted into culture supernatant. The Western blot was detected by using the Flag antibody. (C) Cell surface binding assay showing potential interaction of BTN3A2-S with neurexins. HEK293T cells were transfected with overexpression vector pCMV-Nlgn1-DsRed2, pCMV-Nlgn3-DsRed2, pCMV-Nrxn1β-DsRed2, pCMV-Nrxn1βA4(−)-DsRed2, pCMV-Nrxn3β-DsRed2 or pCMV-DsRed2 (empty vector) alone for 48 h, then were incubated with the culture supernatant of HEK293T overexpressing pCMV-BTN3A2-S-Flag or pCMV-3Tag-8 (green) for 4 h before fixing for visualization. Scale bar, 5 μm. (D) A working model for the unbalanced expression of BTN3A2 in schizophrenia. I: inhibitory synaptic transmission; E: excitatory synaptic transmission.

Discussion

Hitherto, GWAS studies have identified hundreds of loci associated with schizophrenia [7,10]. We now know that the genetic basis and heritability of schizophrenia were largely facilitated by common variations [7,10]. However, as many of these risk variants were located in intergenic regions or other non-coding regions, identifying ‘targetable’ genes in sizable studies has been proven to be difficult, with a few exceptions [40,47,58,59]. Moreover, due to the complexity of schizophrenia and different genetic backgrounds among various ethnic groups, it is even harder to find out the causal genes conferring the risk of schizophrenia. Considering the fact that most of the identified genetic risk loci are located in non-coding regions, these factors may confer the risk of schizophrenia by affecting the expression of nearby gene (s) and/or changing the chromatin state [13,14,60]. Integrative studies with brain gene expression data or other multiple-omics data like DNA methylation, proteins or metabolites can help us parse the result of GWAS studies, as we and others had done recently [[12], [13], [14], [15], [16]]. However, the result of integrative analysis only provided a list of prioritized genes [[12], [13], [14], [15], [16]]. Further follow-up functional experiments were essential to verify the results [13], which constitutes the main work in the post-GWAS era.

In this study, we used a large-scale genotype and expression data from 268 human DLPFC samples from the BrainCloud [26] to conduct eQTL analyses. Then we integrated these eQTL data [26] with CLOZUK + PGC GWAS data [7] results to identify risk genes whose expression changes might confer risk of schizophrenia. We found that BTN3A2 was a potential schizophrenia risk gene, and increased expression of BTN3A2 may contribute to schizophrenia pathogenesis. Further validation using GTEx brain eQTL datasets [32], PsychENCODE brain PFC eQTL datasets [31] and a different integrative analysis method [37] confirmed our findings. Moreover, integrating GWAS data [7] with meQTL data [33] also identified the same gene. These integrative studies from gene expression and methylation suggested that BTN3A2 might affect the risk of schizophrenia, and altered methylation statuses of its nearby sites might play a role. While integrating GWAS and eQTL data is widely accepted as a promising strategy to reveal mechanisms underlying the genetic risk of schizophrenia, few of these studies conducted in recent years have successfully explained how the GWAS risk SNPs affected the expression of nearby genes [13,14,16]. In this study, we combined many newly published datasets like CLOZUK + PGC GWAS data [7] and PsychENCODE brain PFC eQTL data [31] that had not been used in previous schizophrenia integrative studies [13,16,24,61]. Reports from Ma et al. [24], Luo et al. [16] and Yang et al. [13] used the eQTL dataset generated by Myers et al. [62], which was published many years ago. However, in this study, we combined the newly published CLOZUK + PGC GWAS data [7] with RNA-seq generated eQTL dataset [32] using newly developed SMR method [28] and deposited all the integrative analysis results in SZDB (www.szdb.org/integrative1.php) [12] for data sharing. Moreover, we found that the expression and methylation levels of BTN3A2 were significantly changed during the prenatal and postnatal neurodevelopment, with an absolute opposite trend. We further showed that the mRNA levels of BTN3A2 were significantly up-regulated in brain tissues of schizophrenia cases compared with controls. Consistent with the finding in brains, we found that the BTN3A2 gene was also significantly up-regulated in induced pluripotent stem cells derived neurons from schizophrenia patients [41]. Furthermore, overexpression of BTN3A2 in rat hippocampal CA1 pyramidal neurons could decrease both AMPAR and NMDAR-mediated synaptic transmissions, but had no effect on GABA receptor-mediated inhibitory postsynaptic transmission. Finally, cell surface binding assays indicated that this effect might be mediated by the binding of BTN3A2 to neurexins. This study provided a good example of identifying risk gene for schizophrenia via integrative analyses and functional characterization.

Schizophrenia has long been considered as a neurodevelopmental disorder resulted from aberrant developmental processes in fetal, childhood, or adolescent stage [63]. The central nervous system, especially the brain, is vulnerable to the impact of genetic and environmental risk factors for schizophrenia in the prenatal and postnatal developmental periods [63]. Indeed, previous studies have found that many genes that play key roles in neurodevelopment, including brain circuits and synaptic organization, are potential schizophrenia risk genes. For example, ZNF804A is the first putative schizophrenia risk gene identified by GWAS [3,64]. The risk allele of schizophrenia risk SNP rs1344706 was found to decrease ZNF804A cis-expression in the fetal brain during the second trimester, but had no regulatory effects in the first trimester or the adult brain [65]. In utero mouse overexpression experiment verified that ZNF804A was required for normal progenitor proliferation and neuronal migration [66]. In addition, NRGN, another important synapse organizer, could ameliorate ZFP804A-mediated migration defect [66]. In this study, we found that the mRNA expression of BTN3A2 was significantly changed during the prenatal and postnatal neurodevelopment, which might be disturbed by different alleles of GWAS risk SNPs. Importantly, BTN3A2 could interact with presynaptic neurexins and this interaction might lead to a decrease of the excitatory presynaptic release. These lines of evidence further suggested that BTN3A2 might be a schizophrenia susceptibility gene by altering excitatory synaptic function. It should be mentioned that we recently identified independent major depressive disorder (MDD) risk SNPs in the MHC region, which would also tag BTN3A2 as a MDD risk gene [67].

Based on the data from Gentree database (http://gentree.ioz.ac.cn/) [44] and evolutionary analysis, we found that BTN3A2 is specific to some Old World monkeys and human species [45]. Owing to the complexity of cerebral cortex, the human brain developed advanced cognitive, emotional and social capacity. However, the more complex the cerebral cortex is, the more likely that aberrant neural connectivity would appear and result in mental disorders [68,69]. A good case is human-specific NOTCH2NL gene, which contributes to the rapid evolution of a larger human neocortex with the trade-off effect of loss of genomic stability at the 1q21.1 locus, finally confers recurrent neurodevelopmental disorders [70]. Whether this newly evolved BTN3A2 in primates is associated with increased the complexity of the brain, with the cost of higher schizophrenia risk, remains to be answered.

The current study has some limitations. First, the PHEIDI values of SMR [28] integrating CLOZUK + PGC GWAS data [7] with BrainCloud eQTL [26] or CLOZUK + PGC GWAS data [7] with PsychENCODE eQTL data [31] were < 0.01 (Table 1 and Table S3), which suggested that the SMR association was most likely driven by linkage [28,71]. However, the region containing BTN3A2 has a very complex LD pattern (Fig. S4), and it would be hard to identify the causal SNP(s). Second, the risk gene list distilled by integrative analyses was dependent on gene expression profiles of eQTL datasets. Different eQTL datasets might offer different results. For instance, a recent integrative analysis using SMR [28] to integrate CLOZUK + PGC GWAS data [7] with fetal brain eQTL [72] did not identify the association of BTN3A2 with schizophrenia. Evidently, refined expression data would promote the authenticity of the identified risk gene. A recent analysis of brain single cells by single cell RNA-sequencing highlighted a key role of four types of cells (hippocampal CA1 pyramidal cells, striatal medium spiny neurons, neocortical somatosensory pyramidal cells and cortical interneurons) in schizophrenia [49]. A comprehensive integration of expression data at the single-cell level would cast new insights into the etiology of schizophrenia. Third, the electrophysiological analyses of rat brain with overexpression of BTN3A2 might not be a perfect system to show the neurobiological function of this gene as this gene was primate-specific [45]. Tests using non-human primates (NHPs) lacking BTN3A2 and/or overexpression of this gene should be a promising way to approach the final answer, given the recent advancements in genetic modification of NHPs [[73], [74], [75]].

In summary, our integrative analyses and functional assays showed that BTN3A2 is a schizophrenia susceptibility gene. Altered expression of BTN3A2 could affect excitatory and inhibitory synaptic balance through regulating excitatory presynaptic release. Further studies using the schizophrenia related single cell RNA-seq eQTL data and NHP models may help with elucidating the etiology of schizophrenia.