Meta-analysis supports GWAS-implicated link between GRM3 and schizophrenia risk.

Genome-wide association study (GWAS) evidence has identified the metabotropic glutamate receptor 3 (GRM3) gene as a potential harbor for schizophrenia risk variants. However, previous meta-analyses have refuted the association between GRM3 single-nucleotide polymorphisms (SNPs) and schizophrenia risk. To reconcile these conflicting findings, we conducted the largest and most comprehensive meta-analysis of 14 SNPs in GRM3 from a total of 11 318 schizophrenia cases, 13 820 controls and 486 parent-proband trios. We found significant associations for three SNPs (rs2237562: odds ratio (OR)=1.06, 95% confidence interval (CI)=1.02-1.11, P=0.017; rs13242038: OR=0.90, 95% CI=0.85-0.96, P=0.016 and rs917071: OR=0.94, 95% CI=0.91-0.97, P=0.003). Two of these SNPs (rs2237562, rs917071) were in strong-to-moderate linkage disequilibrium with the top GRM3 GWAS significant SNP (rs12704290) reported by the Schizophrenia Working Group of the Psychiatric Genomics Consortium. We also found evidence for population stratification related to rs2237562 in that the 'risk' allele was dependent on the population under study. Our findings support the GWAS-implicated link between GRM3 genetic variation and schizophrenia risk as well as the notion that alleles conferring this risk may be population specific.

Introduction

The Schizophrenia Working Group of the Psychiatric Genomics Consortium (PGC) has identified several significant genome-wide associations between genes involved in glutamatergic neurotransmission and schizophrenia.[1] Among these associations, the metabotropic glutamate receptor 3 (GRM3) contained the second most significant single-nucleotide polymorphism (SNP) (rs12704290, P=3.33 × 10−10).[1] This PGC finding complements a number of previous associations between GRM3 and schizophrenia-related phenotypes, such as prefrontal activation during cognitive tasks,[2] white matter integrity,[3] hippocampal volume in severe obstetric complications,[4] poor cognitive performance (that is, verbal fluency, digit symbol test, perseverative error processing and spatial working memory),[5, 6, 7, 8, 9] antipsychotic response in first-episode schizophrenia[10] and positive symptoms severity in treatment-resistant schizophrenia.[5] However, two meta-analyses of GRM3 SNPs have failed to support an association with schizophrenia risk.[11, 12]

In the first meta-analysis, Albalushi et al.,[11] studied two GRM3 SNPs (rs1468412 and rs2299225) and more recently Yang et al.,[12] examined five GRM3 SNPs (rs1468412, rs274622, rs917071, rs2299225 and rs6465084). Importantly, in both these meta-analyses, the included studies primarily comprised individuals from Asian populations. This is in contrast to the PGC genome-wide association study (GWAS) that predominately included individuals of European descent, suggesting that the association between GRM3 and schizophrenia may be population specific. Furthermore, neither of the previous meta-analyses included family-based genetic association studies, which may have biased their results.

Since the publication of the most recent GRM3 meta-analysis,[12] there have been three additional studies published involving 6027 (2381 schizophrenia and 3646 control) individuals,[2, 6, 13] indicating a reappraisal is warranted. Thus, the aim of this study was to provide an updated and more comprehensive meta-analysis of the association between GRM3 genetic variation and schizophrenia by including both case–control and family studies as well as an exploration of potential population stratification.

Materials and methods

Search strategy

A systematic search of three electronic databases: PubMed, PsycINFO and Medline (Ovid) was performed using combinations of the following keywords: ‘metabotropic glutamate receptor 3’, ‘GRM 3’, ‘mGluR3’ and ‘schizophrenia’. The search was limited to human studies published between January 2002 and December 2016. Two reviewers (SMS and MSM) independently screened the abstracts of all articles identified by the search strategy and then assessed the full-text copies of the eligible articles and extracted the data. In cases where genotype data were unavailable or incomplete, we contacted the corresponding authors and requested the data. A manual search of article bibliographies was completed to identify any non-indexed articles. The SZGene database (www.szgene.org) was also used as a resource for genotype data collection. The ‘preferred reporting items for systematic review and meta-analysis protocols’ (PRISMA-P)[14] was followed in the reporting of this meta-analysis.

Inclusion criteria and data extraction

The following inclusion criteria were used for study selection: (1) examination of GRM3 polymorphisms in schizophrenia, (2) case–control or family-based genetic association study, (3) schizophrenia diagnosis using a standard classification system (for example, DSM-IV, ICD 10 and so on)[15] and (4) provided sufficient genotype or allelic data for calculation of an odds ratio (OR).

From the included studies, the following data were extracted: (a) author(s) and publication year, (b) number of cases and controls or family sample size, (c) country of origin or ethnicity of study participants, (d) diagnostic criteria used, (e) SNP reference sequence number or marker identifier, (f) the publication identification number (for example, PubMed ID) and (g) genotype counts and/or allele counts in cases and controls or family samples (transmission/non-transmission from heterozygous parents to affected offsprings). The extracted data are provided in Supplementary File 1. We also extracted haplotype data when available, but these data were not analyzed because there were no two haplotypes that appeared in three or more selected studies.

Statistical analysis

Data were analyzed using R version 3.3.0 (R Foundation for Statistical Computing, Vienna, Austria). The meta[16] and metafor[17] packages were used to conduct the meta-analyses. The OR with 95% confidence intervals (CIs) was used as the effect size estimator assuming an additive genetic model.

The method proposed by Kazeem and Farrall[18] was used to calculate the effect size for transmission disequilibrium test studies, where the ORs were estimated from the number of transmissions versus non-transmissions of the minor or high-risk allele to schizophrenia cases from heterozygous parents. For case–control studies, ORs were calculated by contrasting the ratio of counts of the minor versus major allele or, when known, high-risk versus low-risk allele in schizophrenia cases versus healthy controls. Two-tailed P-values were reported and were indicated accordingly in the text.

All statistical tests (except for the Q-statistic) were considered statistically significant at P<0.05. Due to differences in study design and sample characteristics, considerable heterogeneity was expected between the studies. Therefore, the pooled ORs were calculated using the random-effects models with the DerSimionian–Laird (DL) estimator,[19] which is based on a normal distribution. The standard error estimates were adjusted using the Hartung–Knapp–Sidik–Jonkman[20, 21] correction, which then calculates the corresponding 95% CIs based on the t-distribution. The Hartung–Knapp–Sidik–Jonkman method generally outperforms the DL approach on type-I error rates when there is heterogeneity, and the number of studies in the meta-analysis is small.[22, 23]

The presence of outliers or influential studies may distort the robustness of the conclusions from the meta-analysis, and as such it is recommended that analysis of potential outlier and influential studies should be an integral part of meta-analysis.[24, 25] Outliers and influential studies were identified according to the recommendations of Viechtbauer and Cheung.[26] Studies with observed effects that are well separated from the rest of the data are considered outliers. Such studies were identified using studentized deleted residuals, with absolute values >1.96 indicative of outliers. An influential study leads to considerable changes to the fitted model and a range of case-deletion diagnostics adapted from linear regression can be used to identify these studies, including the DFFITS, DFBETAS and COVRATIO statistics (see Viechtbauer and Cheung[26] for more information). Potential outliers and influential studies were omitted, and the analyses were then re-run to determine their influence on the pooled effect size. Removal of these studies would substantially reduce the amount of heterogeneity and increase the precision of the pooled effect size, that is, narrowing the confidence interval.

Heterogeneity in effect sizes across studies was tested using the Q-statistic (with P<0.10 indicating significant heterogeneity) and its magnitude was quantified using the I2 statistic, which is an index that describes the proportion of total variation in study effect size estimates that is due to heterogeneity and is independent of the number of studies included in the meta-analysis and the metric of effect sizes.[27] As the Q-statistic has low power when the number of studies is small,[28] 95% prediction intervals were calculated to quantify the extent of heterogeneity in the distribution of effect sizes.[29] The prediction interval is an estimation of the range within which 95% of the true effect sizes are expected to fall.

Moderator analyses for study design (case–control versus transmission disequilibrium test) and ancestry (East Asian versus European) were conducted using mixed-effects meta-analyses. For this method, studies within potential moderator groups were pooled with the random-effects model, while tests for significant differences between the groups were conducted with the fixed-effects model. The Hartung–Knapp–Sidik–Jonkman adjustment was used if there were at least three studies in each group; otherwise, the unadjusted DL method was used.

Publication bias was evaluated by funnel plots and the trim-and-filled procedure,[30] which yields an estimate of the effect size after publication bias has been taken into account.[31] A test for funnel plot asymmetry was only performed if the number of studies was 10 or greater.[32] The test proposed by Harbord et al.[33] was used to quantify the bias captured by the funnel plot and tested whether it was statistically significant.

Results

A total of 20 (17 case–control, 3 family-based) studies representing 11 318 schizophrenia cases, 13 820 controls and 486 parent–proband trios were eligible for further analysis (Supplementary Figure S1 and Table 1). Among the 48 GRM3 SNPs examined in these studies, 14 SNPs were assessed in three or more studies and were subjected to meta-analysis (Figure 1). Significant associations were not found for any of the 14 SNPs when examining all available studies (Supplementary Table S1). However, after removal of an outlier study[34] the rs2237562 C allele was associated with greater schizophrenia risk (k=4, OR=1.06, 95% CI=1.02–1.11, P=0.017) (Figure 2a). In contrast, removal of influential studies[35, 36] showed the C alleles for both rs13242038 (k=3, OR=0.90, 95% CI=0.85–0.96, P=0.016) and rs917071 (k=8, OR=0.94, 95% CI=0.91–0.97, P=0.003) to be associated with lower risk of the disorder (Figures 2b and c).

Table 1

Descriptive characteristics of selected studies in meta-analysis

Author	Year	Case/control	Ancestry	Method	Criteria
Marti et al.[37],a	2002	265/227	German	Case–control	DSM-III-R
Martiet al.[37],a	2002	288/162	German	Case–control	DSM-III-R
Martiet al.[37],a	2002	128	German	Family study	DSM-III-R
Fujiiet al.[38]	2003	100/100	Japanese	Case–control	DSM-IV
Chen et al.[39]	2005	752/752	Chinese	Case–control	DSM-IV
Norton et al.[35]	2005	674/716	Caucasian, UK	Case–control	DSM-IV
Tochigi et al.[36]	2006	402/468	Japanese	Case–control	DSM-IV
Bishop et al.[40]	2007	207/147	Caucasian, US	Case–control	CASH
Schwab et al.[34]	2008	242	German	Family study	DSM-III
Mössner et al.[9]	2008	631/519	German	Case–control	ICD-10
Albalushi et al.[11]	2008	1916/1915	Japanese	Case–control	DSM-IV
Nunokawa et al.[41]	2008	2358/2433	Japanese	Case–control	DSM-IV
Nicodemus et al.[42]	2008	116	Caucasian	Family study	DSM-IV
Betcheva et al.[43]	2009	255/556	Bulgarian	Case–control	DSM-IV
Jönsson et al.[44]	2009	582/1473	Scandinavian	Case–control	DSM-IV
Jia et al.[45]	2014	433/619	Chinese	Case–control	DSM-IV
Mounce et al.[3]	2014	74/87	Caucasian, US	Case–control	DSM-IV
O'Brien et al.[13]	2014	1235/1309	Caucasian, UK	Case–control	DSM-IV
Chang et al.[6]	2015	1115/2289	Chinese	Case–control	DSM-IV
Kinoshita et al.[2]	2015	31/48	Japanese	Case–control	DSM-IV

These three studies were published in the same article.[37]

Figure 1

Schematic diagram of the GRM3 structure with locations of SNPs is examined in this meta-analysis. Single-nucleotide polymorphisms (SNPs) significantly associated with schizophrenia in this meta-analysis are underlined. *This SNP was excluded in this meta-analysis due to fewer than three studies but was the most significant GRM3 SNP in the PGC GWAS.[1]

Figure 2

Forest plots of (a) rs2237562, (b) rs13242038 and (c) rs917071 association with schizophrenia before and after excluding outlier/influential studies. CI, confidence interval; OR, odds ratio; RE model, random-effect model.

Low to moderate heterogeneity (I2=0.0–34.6%) was present across the studies examining these three significant SNPs, although removal of outlier/influential studies reduced heterogeneity to 0% for all three SNPs (Supplementary Table S1). No evidence of moderation by study design or ancestry was detected (Supplementary Table S2 and S3).

Publication bias was evaluated for rs6465084 and rs1468412 as these were the only two SNPs with ten or more studies. The regression test for funnel plot asymmetry was statistically significant for studies examining rs1468412 (t[10]=2.35, P=0.040) (Supplementary Figure S2). After adjustment for publication bias with the trim-and-fill procedure, the odds ratio was decreased to 0.96 (95% CI=0.86, 1.07, P=0.426; the number of filled studies was 2. Given that there was a subtle reduction in the pooled odds ratio, the impact of publication bias was likely modest. For rs6465084 the regression test for funnel plot asymmetry was not statistically significant (t[8]=0.58, P=0.578) (Supplementary Figure S3).

Discussion

We found three SNPs (rs2237562, rs13242038, rs917071) in GRM3 that were significantly associated with schizophrenia. Our results differ from the two previous meta-analyses[35, 36] that did not detect an association between GRM3 genetic variation and schizophrenia. However, our positive results were only apparent after removal of outlier or influential studies, a method recommended[24, 25] but not employed in the previous two meta-analyses.[11, 12]

Importantly, our meta-analysis supports finding from the PGC GWAS, which identified GRM3 as one of 108 loci associated with schizophrenia.[1] In fact, the top PGC GWAS SNP in GRM3 (rs12704290, P=1.0 × 10−10) is in moderate to strong linkage disequilibrium (LD) with two (rs2237562: D′=1.00, r2=0.485; rs917071: D′=0.64, r2=0.199) of the three significant SNPs we detected in the current meta-analysis. The association between rs2237562 and schizophrenia is particularly interesting for two reasons. First, among the 14 SNPs eligible for meta-analysis, the rs2237562 SNP had the strongest LD with the top GRM3 GWAS hit (rs12704290), which has only been examined in one case–control study[6] to date. Second, in the PGC GWAS the T allele of rs2237562 was associated with schizophrenia risk (OR=1.04, P=8.3 × 10−5) while in our meta-analysis the C allele was associated with risk. One potential explanation for this discrepancy is population stratification. The PGC GWAS was predominately focused on people of European descent with people of Asian descent represented in only 3.5% (Japanese: 0.6%, Singapore: 1.2% and Chinese: 1.7%) of the studied population. Conversely, the four studies contributing to our significant meta-analysis finding were exclusively conducted in the East Asian population.[6, 36, 39, 45] Notably, the two alleles for rs2237562 (T/C) occurred within similar haplotype blocks and at comparable frequencies in the European and in East Asian population (Supplementary Figures S4 and S5).

In fact, inclusion of the one outlier study,[34] which was conducted in a German population and concurred with the PGC findings, eliminated the pooled effect detected among the more ancestral homogenous East Asian studies. Thus, the ‘risk’ allele for rs2237562 may be ancestry dependent, with potential implications for calculating polygenic risk scores that include this SNP within non-European populations.

It is also noteworthy that current evidence supporting an association between GRM3 genetic variation and schizophrenia susceptibility as well as cognitive and neuroimaging indices have been centered on SNPs in exon 3 and its adjacent introns.[46] However, among the five SNPs we examined in this region, only the intronic rs2237562 SNP (discussed above) was supported. Yet, upstream of this region in the intron between exons I and II, we showed the C alleles for two SNPs (rs13242038 and rs917071) were associated with schizophrenia risk. Conversely, the alternative (T) allele for both SNPs was associated with lower schizophrenia risk (rs13242038: OR=0.95, P=2.7 × 10−4; rs917071: OR=0.95, P=0.056) in the PGC GWAS, suggesting again a potential population stratification effect. However, the studies contributing to the meta-analyses of these two SNPs were not ancestrally homogenous (rs13242038: 68% East Asian and 32% Caucasian; rs917071: 49% East Asian and 51% Caucasian) and moderation analysis did not detect an ancestry effect. Furthermore, the removed influential studies by Norton et al.[35] for rs13242038 and Tochigi et al.[36] for rs917071 were conducted in Caucasian and Asian populations, respectively. Thus, the discordance between our meta-analytical findings and those of the PGC is likely not a result of population structure. Previous research supports the PGC findings, in part, in that individuals with schizophrenia who had the rs13242038 T allele and a history of severe obstetric complications (OCs) had larger hippocampi. However, similar support for rs917071 is absent in that previous neuroimaging[4] and cognitive[6, 7] studies in schizophrenia have been negative to date.

The mechanisms underpinning the genetic associations between the three GRM3 SNPs and schizophrenia are still unknown. Gene expression, more specifically gene splicing, is postulated to mediate the pathophysiological effects of allelic variation.[46, 47] Four alternative splice variants have been reported in the human brain: full-length GRM3 (2.8 kb), GRM3 with exon 2 deleted (GRM2Δ3, 2.2 kb), GRM3 with exon 4 deleted (GRM3Δ4, 1.8 kb) and GRM3 with exons 2 and 3 deleted (GRM3Δ2Δ3, 1.4 kb).[48, 49] A previous study[50] found that an exon 3 SNP (rs2228595) predicted increased transcript expression of the GRM3Δ4 splice variant that encodes a truncated receptor in the dorsolateral prefrontal cortex of schizophrenia patients, suggesting GRM3 SNPs that are in LD with rs2228595 may also affect GRM3 splicing. Interestingly, the rs2228595 is in strong LD with the PGC GWAS SNP (rs12704290: D′=1.0, r2=0.011) and our meta-analysis significant SNPs (rs2237562: D′=1.0, r2=0.170; rs917071: D′=0.757, r2=0.098), but we found no association of rs2228595 with schizophrenia in the current meta-analysis (k=5, OR=1.06, 95% CI=0.73–1.54, P=0.706). Nevertheless, there is a preliminary indication that GRM3 genetic variation may influence gene splicing, but additional research interrogating how this splicing may contribute to the pathophysiology of schizophrenia is warranted.

In conclusion, we have identified novel meta-analytically derived associations between polymorphisms within the fourth (rs2237562) and second (rs13242038, rs917071) introns of GRM3 with schizophrenia. Two of these SNPs (rs2237562, rs917071) were in strong-to-moderate LD with the top PGC SNP identified in the largest schizophrenia GWAS to date.[1] We further showed that the rs2237562 allele associated with ‘risk’ for schizophrenia is likely ancestry dependant. As such, our results support the GWAS-implicated link between GRM3 genetic variation and schizophrenia risk as well as support the notion that alleles conferring this risk may be population specific. These findings provide further justification for the future study of the mechanism(s) by which GRM3 genetic variation contributes to schizophrenia risk.