Advertisement

Rare-variant pathogenicity triage and inclusion of synonymous variants improves analysis of disease associations of orphan G protein–coupled receptors

Open AccessPublished:October 18, 2019DOI:https://doi.org/10.1074/jbc.RA119.009253
      The pace of deorphanization of G protein–coupled receptors (GPCRs) has slowed, and new approaches are required. Small molecule targeting of orphan GPCRs can potentially be of clinical benefit even if the endogenous receptor ligand has not been identified. Many GPCRs lack common variants that lead to reproducible genome-wide disease associations, and rare-variant approaches have emerged as a viable alternative to identify disease associations for such genes. Therefore, our goal was to prioritize orphan GPCRs by determining their associations with human diseases in a large clinical population. We used sequence kernel association tests to assess the disease associations of 85 orphan or understudied GPCRs in an unselected cohort of 51,289 individuals. Using rare loss-of-function variants, missense variants predicted to be pathogenic or likely pathogenic, and a subset of rare synonymous variants that cause large changes in local codon bias as independent data sets, we found strong, phenome-wide disease associations shared by two or more variant categories for 39% of the GPCRs. To validate the bioinformatics and sequence kernel association test analyses, we functionally characterized rare missense and synonymous variants of GPR39, a family A GPCR, revealing altered expression or Zn2+-mediated signaling for members of both variant classes. These results support the utility of rare variant analyses for identifying disease associations for GPCRs that lack impactful common variants. We highlight the importance of rare synonymous variants in human physiology and argue for their routine inclusion in any comprehensive analysis of genomic variants as potential causes of disease.

      Introduction

      The superfamily of G protein–coupled receptors (GPCRs)
      The abbreviations used are: GPCR
      G protein–coupled receptor
      GWAS
      genome-wide association study
      SKAT
      sequence kernel association test
      LOF
      loss of function
      MISS
      missense
      EHR
      electronic health record
      MAF
      mean allele frequency
      pLP
      predicted likely pathogenic
      pP
      predicted pathogenic
      PheWAS
      phenome-wide association studies
      IP3
      inositol trisphosphate
      IP1
      inositol monophosphate
      CI
      confidence interval
      pB
      predicted to be benign
      WES
      whole exome sequences
      ASCVD
      atherosclerotic cardiovascular disease.
      translates extracellular signals from light, metabolites, and hormones into cellular changes, which makes them the targets of a significant fraction of drugs currently on the market (
      • Hauser A.S.
      • Attwood M.M.
      • Rask-Andersen M.
      • Schiöth H.B.
      • Gloriam D.E.
      Trends in GPCR drug discovery: new agents, targets and indications.
      ). Genome-wide association studies (GWASs) on common variants in GPCRs have begun to identify their contributions to various disease processes (
      • Roth B.L.
      • Kroeze W.K.
      Integrated approaches for genome-wide interrogation of the druggable non-olfactory G protein–coupled receptor family.
      ,
      • Kovacs P.
      • Schönberg T.
      The relevance of genomic signatures at adhesion GPCR loci in humans.
      ). However, many GPCRs lack functionally important common variants, and alternate strategies are needed to understand their roles in human physiology. Recently, sequence kernel association tests (SKATs) have been applied to rare variants to assess disease associations. Initial approaches binned all rare variants in a genomic region or gene (
      • Wu M.C.
      • Lee S.
      • Cai T.
      • Li Y.
      • Boehnke M.
      • Lin X.
      Rare-variant association testing for sequencing data with the sequence kernel association test.
      ,
      • Lee S.
      • Emond M.J.
      • Bamshad M.J.
      • Barnes K.C.
      • Rieder M.J.
      • Nickerson D.A.
      • NHLBI GO Exome Sequencing Project–ESP Lung Project Team
      • Christiani D.C.
      • Wurfel M.M.
      • Lin X.
      Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies.
      ), but current methods group the rare variants most likely to contribute to disease associations or aggregate variants based on domain or family structures (
      • Richardson T.G.
      • Shihab H.A.
      • Rivas M.A.
      • McCarthy M.I.
      • Campbell C.
      • Timpson N.J.
      • Gaunt T.R.
      A protein domain and family based approach to rare variant association analysis.
      ,
      • Friedrichs S.
      • Malzahn D.
      • Pugh E.W.
      • Almeida M.
      • Liu X.Q.
      • Bailey J.N.
      Filtering genetic variants and placing informative priors based on putative biological function.
      ).
      Large-scale studies of orphan or understudied GPCRs to characterize natural genetic variation in the human population can provide insights into the biological function and/or potential causal contributions to disease processes (
      • Sriram K.
      • Insel P.A.
      G protein–coupled receptors as targets for approved drugs: how many targets and how many drugs?.
      ). The DiscovEHR collaboration represents a tremendous resource (
      • Dewey F.E.
      • Murray M.F.
      • Overton J.D.
      • Habegger L.
      • Leader J.B.
      • Fetterolf S.N.
      • O'Dushlaine C.
      • Van Hout C.V.
      • Staples J.
      • Gonzaga-Jauregui C.
      • Metpally R.
      • Pendergrass S.A.
      • Giovanni M.A.
      • Kirchner H.L.
      • Balasubramanian S.
      • et al.
      Distribution and clinical impact of functional variants in 50,726 whole-exome sequences from the DiscovEHR study.
      ) because whole exome sequences can be linked to the electronic health record (EHR) in a deidentified manner. Here we used whole exome sequences and clinical information from 51,289 individuals in the DiscovEHR cohort. We identified sequence variants (common variants, mean allele frequency (MAF) ≥ 1%; and rare variants, MAF < 1%) for 85 orphan or understudied GPCRs not amenable to GWAS studies. Rare variants were binned into functional classes: loss-of-function (LOF, truncation and frameshift) variants; missense variants with predicted pathogenicity; or rare synonymous variants with altered codon bias. We performed independent SKAT analyses on the various functional classes to determine their disease associations. Remarkably, for those GPCRs with sufficient numbers of individuals for statistical analyses, we found that the top disease associations were common to all available variant classes in the GPCR.
      For initial in vitro assessment of the validity of the disease associations, we focused on GPR39. Among its particular advantages, GPR39 had sufficient numbers of variants in all three functional classes to be amenable to rare variant approaches. Second, GPR39 is a small family A GPCR, allowing rapid generation of mutants. Third, GPR39 is activated by extracellular Zn2+ and is coupled to inositol phosphate production (
      • Asraf H.
      • Salomon S.
      • Nevo A.
      • Sekler I.
      • Mayer D.
      • Hershfinkel M.
      The ZnR/GPR39 interacts with the CaSR to enhance signaling in prostate and salivary epithelia.
      ,
      • Cohen L.
      • Asraf H.
      • Sekler I.
      • Hershfinkel M.
      Extracellular pH regulates zinc signaling via an Asp reside of the zinc-sensing receptor (ZnR/GPR39).
      ), permitting straightforward functional analyses. Finally, we have only a rudimentary understanding of the role(s) of GPR39 in human physiology, despite its broad expression and multiple cellular functions (
      • Sunuwar L.
      • Gilad D.
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ,
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ); thus, results may provide insights into the contribution(s) of GPR39 variation to disease. Our results demonstrate the validity of combined computational and functional approaches to provide important insights into the clinical contributions of orphan and understudied GPCRs, and focus attention on the importance of rare synonymous variants for identifying disease associations. The overall strategy can be readily applied to other classes of genes without clinically important common variants.

      Results

      Rare variant-based disease associations of orphan or understudied GPCRs

      We identified rare variants in whole exome sequences from 51,289 individuals. To predict their functional impact, variants were annotated and sorted into three classes in order of predicted severity, i.e. LOF (premature stop codon, loss of a start or stop codon, or disruption of a canonical splice site), missense variants (predicted to be likely pathogenic (pLP) or pathogenic (pP) by RMPath scoring, derived from application of multiple pathogenicity algorithms (see “Experimental procedures”)), and synonymous variants with significantly different codon frequency relative to the reference codon (termed SYN_ΔCB, i.e. synonymous variants with altered codon bias).
      Determining disease associations for genes having only low-frequency variants requires binning of variants to increase statistical power. Binning methods have focused on gene or genetic region and have recently been expanded to incorporate regulatory regions and/or pathways by incorporating biological information from curated knowledge databases (
      • Basile A.O.
      • Byrska-Bishop M.
      • Wallace J.
      • Frase A.T.
      • Ritchie M.D.
      Novel features and enhancements in BioBin, a tool for biologically inspired binning and association analysis of rare variants.
      ). For this study, we were specifically looking for clinically impactful rare coding variants that could be validated by functional studies of the relevant GPCRs. We binned variants according to the functional classes described above and used SKAT, focusing on the disease associations that were common to more than one class of rare variant. Fig. 1 illustrates the analysis pipeline for the entire set of 85 GPCRs. Of the 72 GPCRs with sufficient variants to allow at least 2 variant classes to be analyzed by SKAT, 28 yielded significant disease associations. Because this was an exploratory analysis, we used a significance cutoff of p < 0.05, i.e. p values < 6.2E-07 after Bonferonni correction for multiple testing. Table 1 shows the top disease associations for the 13 GPCRs that had sufficient variants in the LOF class. Supporting validity, associations between GPR37L1 and epilepsy (
      • Giddens M.M.
      • Wong J.C.
      • Schroeder J.P.
      • Farrow E.G.
      • Smith B.M.
      • Owino S.
      • Soden S.E.
      • Meyer R.C.
      • Saunders C.
      • LePichon J.B.
      • Weinshenker D.
      • Escayg A.
      • Hall R.A.
      GPR37L1 modulates seizure susceptibility: evidence from mouse studies and analyses of a human GPR37L1 variant.
      ) and LGR5 and alopecia (
      • Hoeck J.D.
      • Biehs B.
      • Kurtova A.V.
      • Kljavin N.M.
      • de Sousa E.
      • Melo F.
      • Alicke B.
      • Koeppen H.
      • Modrusan Z.
      • Piskol R.
      • de Sauvage F.J.
      Stem cell plasticity enables hair regeneration following LGR5+ cell loss.
      ,
      • Michel L.
      • Reygagne P.
      • Benech P.
      • Jean-Louis F.
      • Scalvino S.
      • Ly Ka So S.
      • Hamidou Z.
      • Bianovici S.
      • Pouch J.
      • Ducos B.
      • Bonnet M.
      • Bensussan A.
      • Patatian A.
      • Lati E.
      • Wdzieczak-Bakala J.
      • et al.
      Study of gene expression alteration in male androgenic alopecia: evidence of predominant molecular signaling pathways.
      ) have been previously identified by other means. Attesting to the discovery potential, other disease associations were novel, including GPR161 with hyperhidrosis, LGR4 with Sicca syndrome, GPR153 with hyperpotassemia, GPR84 with anemia, and GPR39 with peripheral nerve damage, among others. Notably, the phenome-wide disease associations for the predicted LOF class showed congruence across the other variant classes for the subset of GPCRs having sufficient missense or SYN_ΔCB variants for independent analyses. Although truncation and frameshift variants are easily identified and likely to have the greatest functional impact, missense variants classed as pLP + pP can also have a significant impact on GPCR function. For the subset of GPCRs without sufficient LOF variants, we found significant disease associations for the missense and SYN_ΔCB variants (Table 1). Some of the associations validate previous studies, including GPR183 and disorders of liver (
      • Guillemot-Legris O.
      • Mutemberezi V.
      • Cani P.D.
      • Muccioli G.G.
      Obesity is associated with changes in oxysterol metabolism and levels in mice liver, hypothalamus, adipose tissue and plasma.
      ) and GPR85 with acute myocardial infarction (
      • Li Y.
      • Dong J.
      • Xuan Y.-H.
      • Liu S.-S.
      • Luo J.-Y.
      • Xiao Y.-J.
      • Tian Z.P.
      • Sun Z.J.
      Identification of microRNA expression in a rat model of post-infarction heart failure.
      ). Other associations were novel and represent potential targets for future study, including GPR132 and disturbances in sulfur-bearing amino acid metabolism, GPR176 with asthma, GPR12 with epilepsy, and GPRC5D with renal osteodystrophy. Of note, the SYN_ΔCB variant class had significant disease associations in common with LOF and/or missense variants and most consistently had sufficient variants to be included in the analysis. Parallel SKAT analyses of the three distinct classes of rare variants produced two important results. First, independent analysis of distinct classes of variants yielded concordant disease associations, increasing confidence in the results. Second, the most consistent source of phenome-wide disease associations was SYN_ΔCB. Altogether, these results provide a strong rationale for including all three functional variant classes in disease association analyses.
      Figure thumbnail gr1
      Figure 1Flow diagram for orphan/understudied GPCR SKAT analysis. Rare variants in 85 GPCRs were identified in the 51,289 DiscovEHR WES data set and sorted into LOF, missense, and synonymous categories. Missense variants were further triaged to identify those predicted to be pLP or pP by RMPath scoring and the subset of synonymous variants with large changes in codon usage bias (SYN_ΔCB; see “Experimental procedures”). GPCRs with at least 2 categories of variant having ≥7 distinct members were further analyzed (72 of 85 GPCRs). SKAT analysis was run independently on 123 sets of variants against 656 disease codes having ≥200 individuals with ≥3 independent encounters in the EHR.
      Table 1Orphan GPCRs with congruent disease associations for multiple classes of rare variants
      GPCRICD9DescriptionLOFMISSSYN_ΔCB
      LGR4710.2Sicca syndrome1.58E-311.28E-132.31E-14
      GPR150372.14Other chronic allergic conjunctivitis2.59E-311.45E-325.48E-07
      GPR153276.7Hyperpotassemia2.32E-308.73E-161.14E-22
      GPR149374.3Ptosis of eyelid5.89E-304.62E-081.53E-07
      GPR161780.8Generalized hyperhidrosis3.78E-26×9.03E-27
      GPR37L1345.1Generalized convulsive epilepsy, without intractable epilepsy2.04E-181.92E-194.95E-20
      GPR84280.8Other specified iron deficiency anemias2.31E-152.99E-082.07E-15
      GPR39354.2Lesion of ulnar nerve8.04E-151.62E-098.69E-15
      LGR6583.81Nephritis and nephropathy, not specified as acute or chronic4.52E-11×1.31E-19
      LGR5704Alopecia NOS5.73E-13×1.36E-12
      GPR1287.5Thrombocytopenia NOS1.5E-091.26E-071.5E-09
      GPR18583.81Nephritis and nephropathy, not specified as acute or chronic2.61E-081.01E-083.77E-09
      GPR179270.4Disturbances of sulfur-bearing amino acid metabolism1.38E-07×1.05E-07
      GPR176493.91Asthma, unspecified type, with status asthmaticus×1.35E-336.03E-26
      GPR132216.3Disturbances of sulfur-bearing amino acid metabolism×4.13E-309.61E-24
      GPR85410.1Acute myocardial infarction of other anterior wall×2.06E-261.33E-35
      GPR183573.8Other specified disorders of liver×1.34E-191.5E-19
      GPR12345.1Generalized convulsive epilepsy, without intractable epilepsy×5.64E-145.48E-14
      GPR135350.1Trigeminal neuralgia×1.91E-121.99E-12
      GPR26372.14Other chronic allergic conjunctivitis×4.56E-124.35E-12
      MRGPRE401.1Benign essential hypertension×2.17E-108.61E-08
      GPR15455.6Hemorrhoids NOS×1.61E-101.97E-10
      GPR83250.8Diabetes mellitus type II (non–insulin-dependent)×1.29E-095.71E-12
      GPRC5D588Renal osteodystrophy×7.5E-092.0E-20
      GPCR5B381.81Dysfunction of Eustachian tube×1.4E-071.23E-07
      GPR141786.07Wheezing×3.17E-092.76E-10
      GPR6214.1Lipoma of other skin and subcutaneous tissue×1.82E-072.51E-09
      GPR3141.02Group B Streptococcus infection×2.04E-072.04E-07

      GPR39 common and rare variants in the 51,289 WES DiscovEHR cohort

      Disease associations provide a rich source of hypotheses regarding the potential contributions of GPCRs to human physiology. The first step in defining a causal role for these rare variants, however, is functional validation of the SKAT results. Here we focused in particular on rare synonymous variants with altered codon frequency (SYN_ΔCB), whose functional consequences are not well-understood. Many of the GPCRs have no known agonist, and their dominant signaling pathways have not been characterized (Table 1). We therefore chose to focus on GPR39, the Zn2+ receptor, which has been characterized in cellular and knockout mouse models and for which both agonist and dominant signaling pathways are known (
      • Sunuwar L.
      • Gilad D.
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ,
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ).
      All common and rare variants in GPR39 were identified in 51,289 individuals in the DiscovEHR cohort (
      • Dewey F.E.
      • Murray M.F.
      • Overton J.D.
      • Habegger L.
      • Leader J.B.
      • Fetterolf S.N.
      • O'Dushlaine C.
      • Van Hout C.V.
      • Staples J.
      • Gonzaga-Jauregui C.
      • Metpally R.
      • Pendergrass S.A.
      • Giovanni M.A.
      • Kirchner H.L.
      • Balasubramanian S.
      • et al.
      Distribution and clinical impact of functional variants in 50,726 whole-exome sequences from the DiscovEHR study.
      ) (Fig. 2). None of the GPR39 variants were found in the ClinVar database. Accordingly, all LOF, missense, and synonymous variants were analyzed and classified as described in Table 1. The common variants were used for Phenome-wide Association Studies (PheWAS) and included a missense variant (2:133174764, pA50V) predicted to be benign, a synonymous variant (2:133174999; pT128T), which changed a common to the rarest codon (acA (28%) to acG (12%), i.e. a 2.33-fold decrease in frequency) and a noncoding variant (2:133402607). Notably, no phenome-wide disease associations were detected (Fig. 3A).
      Figure thumbnail gr2
      Figure 2Pathogenicity triage of 173 coding variants in GPR39 identified in 51,289 WES. Common variants with MAF ≥ 1% were categorized as likely benign (MAF between 1 and 5%) and benign (MAF ≥ 5%). All variants having MAF of <1% were sorted to synonymous, LOF (frameshift or premature stop codon), and missense. Missense variants were thus variants of unknown significance (VUS) and subjected to bioinformatics pathogenicity triage using RMPath scoring (“Experimental procedures”) and classified as pB, pLB, pLP, or pP and the subset of SYN variants having a significant change in codon usage bias (SYN_ΔCB).
      Figure thumbnail gr3
      Figure 3Manhattan plots of disease association analyses of GPR39 variants. A, PheWAS analyses (common SNP versus disease phenotypes in the EHR) were run for the indicated GPR39 common variants. No association results exceeded the minimal p value threshold of 0.05, uncorrected for multiple testing. B, rare LOF plus MISS variants were combined (39 total variants), and disease associations were determined by SKAT. Top phenome-wide associations are labeled. C, disease associations of rare SYN_ΔCB variants (21 total variants) were determined by SKAT; top phenome-wide associations are labeled. Full details of associations common to both variant classes (B and C) are presented in . For both B and C, the significance threshold was set at 0.05 (adjusted p value 3.8E-05, 2 sets, 656 codes) and is indicated by a dotted line.
      The rare variants used for SKAT, i.e. LOF, missense (pLP + pP), and SYN_ΔCB variants, were distributed throughout the coding region of GPR39 (Fig. 4). Examination of the 28 missense variants that were categorized as pLP and pP confirmed that the RMPath pipeline identified variants which are likely to impact GPR39 function (Fig. 4). The three pP variants (S78L, C108R, and R133G) included a cysteine that participates in a disulfide bond at the extracellular face of GPR39 (
      • Storjohann L.
      • Holst B.
      • Schwartz T.W.
      A second disulfide bridge from the N-terminal domain to extracellular loop 2 dampens receptor activity in GPR39.
      ). The pLP variants included one that contributes to the Zn2+-binding site (H19R) (
      • Storjohann L.
      • Holst B.
      • Schwartz T.W.
      Molecular mechanism of Zn2+ agonism in the extracellular domain of GPR39.
      ), two variants in a disulfide bond cysteine (C191G and C191W) (
      • Storjohann L.
      • Holst B.
      • Schwartz T.W.
      A second disulfide bridge from the N-terminal domain to extracellular loop 2 dampens receptor activity in GPR39.
      ), a residue in the highly conserved NPXXY motif common to family A GPCRs (N340S) (
      • Yuan S.
      • Filipek S.
      • Palczewski K.
      • Vogel H.
      Activation of G-protein-coupled receptors correlates with the formation of a continuous internal water pathway.
      ), and a variant distal to a putative palmitoylation site in the proximal C terminus (R352Q). Other residues in the pLP category introduce charges or proline residues into transmembrane helices or extracellular loops (S49N, M87I, T122M, A167S, P179T, R193C, S222Y, T291I, A293P, R302Q, A307V, P310A, Y314R, R320Y, S336R, R352Q, R362C, H367R, R386C, and R390C). We also identified the subset of 21 variants that introduced at least a 2-fold change in site-specific codon frequency (SYN_ΔCB; denoted by asterisks in Fig. 4).
      Figure thumbnail gr4
      Figure 4GPR39 topology indicating rare variants identified in 51,289 WES. The GPR39 sequence was color-coded to indicate reference codon usage frequency. For all amino acids having multiple potential codons, red indicates the most common codon, with progression from common to rare colored as gradations to light pink (rare). Asterisks indicate the locations of rare synonymous variants exhibiting the largest changes in codon bias, either from common to rare or rare to common. indicates rare LOF variants (truncations or frameshifts), and @ indicates locations of rare missense variants combined for SKAT analysis.

      Disease associations of GPR39 rare variant classes

      The initial SKAT screen for GPR39 disease associations compared top associations for the three independent classes of rare variants. To more fully explore the results for GPR39 and increase power to identify all significant disease associations, we combined all individuals heterozygous for rare LOF, pLP, and pP variants (termed LOF/MISS). Analysis was restricted to only those ICD9 codes with ≥200 individuals independent of genotype (
      • Carey D.J.
      • Fetterolf S.N.
      • Davis F.D.
      • Faucett W.A.
      • Kirchner H.L.
      • Mirshahi U.
      • Murray M.F.
      • Smelser D.T.
      • Gerhard G.S.
      • Ledbetter D.H.
      The Geisinger MyCode community health initiative: an electronic health record-linked biobank for precision medicine research.
      ), and the results are plotted against ICD9 codes in Fig. 3B. SKAT analysis was also run on rare SYN_ΔCB variants, plotted in Fig. 3C. Both the LOF/MISS and SYN_ΔCB analyses had a significant number of phenome-wide associations in common. All results for both variant classes that were above the significance cutoff for multiple testing (dotted lines in Fig. 3, B and C) are listed in Table 2. Top associations for GPR39 included lesion of the ulnar nerve and benign prostate hyperplasia with lower urinary tract obstruction. Note that individuals who had both rare LOF/MISS and SYN_ΔCB variants were excluded from analyses, because the sequence data are not phased, and we therefore cannot distinguish between 1 allele having both variants versus a variant on each allele. Also of note in Table 2, some significant disease associations were specific to either the LOF/MISS or SYN_ΔCB variant classes. To assess the importance of rare variant pathogenicity triage for binning prior to SKAT analysis, we combined all rare (MAF < 1%) variants (LOF + missense + synonymous = 168 variants) and ran SKAT. The only significant disease association was upper quadrant abdominal pain (ICD9 789.02; p value 8.98E-06; significance cutoff at p < 7.6E-05). We also ran SKAT on those rare missense variants predicted to be benign (Fig. 2, pB; 61 variants); no statistically significant disease associations were obtained (data not shown). Likewise, we ran SKAT analysis on all rare synonymous variants and found top associations similar to those obtained for SYN_ΔCB but with less significant p values, suggesting dilution of impactful variants, i.e. lesion of ulnar nerve (p value 1.08E-07), and benign prostate hyperplasia with urinary tract obstruction and other lower urinary symptoms (p value 3.79E-05) (compare Table 2, Table 3). We conclude that pathogenicity triage to refine rare variant binning was crucial to identifying disease associations for the LOF and missense variants, and binning of the subset of rare synonymous variants having large changes in local codon frequency bias strengthened disease associations. Further, refining the rare variant input into the SKAT analyses yielded disease associations common to multiple independent classes of variants, strengthening confidence in the associations.
      Table 2SKAT results for rare GPR39 variants in 51,289 DiscovEHR cohort
      CodeDescriptionLOF + pLP + pPSYN_ΔCB
      354.2Lesion of ulnar nerve6.2E-088.69E-15
      995.3Allergy NOS6.67E-067.22E-14
      600.21Benign prostate hyperplasia w urinary obstruction1.24E-072.57E-09
      362.56Macular puckering of retina×7.65E-07
      786.6Swelling, mass, or lump in chest3.83E-052.4E-07
      704.8Other specified diseases of hair and hair follicles×6.72E-07
      696Psoriatic arthropathy3.24E-131.03E-06
      727.05Other tenosynovitis of hand or wrist×1.78E-06
      372.3Conjunctivitis NOS1.68E-052.75E-06
      998.59Other postoperative infection×5.0E-06
      401.1Benign essential hypertension×9.55E-06
      70.54Chronic hepatitis C without mention of hepatic coma×1.09E-05
      655.83Other known or suspected fetal abnormality9.71E-061.41E-05
      591Hydronephrosis×2.06E-05
      110.4Dermatophytosis of foot1.62E-06×
      726.2Other affections of shoulder region NOS4.54-E06×
      429.2Cardiovascular disease (ASCVD)2.26E-05×
      707.19Ulcer of lower limb2.46E-05×
      Table 3SKAT analysis of all rare synonymous GPR39 variants in 51,289 DiscovEHR cohort
      CodeDescriptionp value
      696Psoriatic arthropathy3.87E-11
      354.2Lesion of ulnar nerve1.08E-07
      110.4Dermatophytosis of foot8.34-E07
      726.2Other affectations of shoulder region, NOS2.24E-06
      995.3Allergy, NOS3.51E-06
      655.83Other known or suspected fetal abnormality4.38E-06
      707.19Ulcer of lower limb1.06E-05
      786.6Swelling, mass, or lump in chest2.06E-05
      429.2Cardiovascular disease (ASCVD)3.03E-05
      649.13Obesity complicating pregnancy, childbirth, antepartum3.58E-05
      600.21Benign prostate hyperplasia with lower urinary tract symptoms3.79E-05
      440.9Atherosclerosis NOS6.31E-05
      649Tobacco use disorder complicating pregnancy, childbirth6.8E-05
      368.2Diplopia7.46E-05

      Validation of the functional impact of rare GPR39 variants

      Point mutations can cause changes in GPR39 function by altering protein expression and/or by altering essential features of receptor activation. Therefore, to explore the functional basis for the observed disease associations, we generated the rare missense and SYN_ΔCB variants in the reference human GPR39 background. The cDNA construct included an N-terminal 3× FLAG epitope followed by a minimal bungarotoxin-binding site, for ease of expression analysis and subcellular localization, respectively. Expression in HEK293 cells of equivalent amounts of cDNA and subsequent Western blotting of lysates showed significant differences in net expression levels of the rare missense and SYN_ΔCB variants. Fourteen missense variants had reduced and one had significantly increased expression relative to WT GPR39 (Fig. 5, A and C). Rare missense variants with altered expression were distributed throughout GPR39 structural domains as color-coded in the quantitation of multiple blots (Fig. 5C). Likewise, SYN_ΔCB variant expression was significantly greater than WT GPR39 for three variants, whereas expression of fourteen variants was significantly less that WT (Fig. 5B), as quantified in Fig. 5D. Surprisingly, despite the change from rare to common codon frequency for the majority of SYN_ΔCB variants, only four showed increased expression relative to WT. Four variants went from common to rare codons, denoted by asterisks on the x axis in Fig. 5D. Contrary to expectations, Gln-265 was expressed at significantly higher levels than WT (∼150%), despite the large decrease in codon frequency. Overall, results argue that factors in addition to site-specific codon usage bias determine net expression levels for rare synonymous variants.
      Figure thumbnail gr5
      Figure 5Protein expression of rare missense and synonymous GPR39 variants. A, representative Western blots of rare missense variants identified in the 51,289 WES cohort. Western blots of HEK293 cells lysates expressing WT or rare variant missense cDNAs were run as described under “Experimental procedures.” Each blot had HEK and WT for background subtraction and normalization, respectively. The blots were probed for GAPDH to confirm equal loading. B, representative blots for HEK293 cells expressing WT and rare SYN_ΔCB variants. Expression levels varied and to ensure quantitation relative to WT, blots were segregated for low (upper) and high (lower) expressing variants. The blots were cut (dotted line) to allow simultaneous probing for GPR39 and GAPDH. C, quantitation of expression levels of GPR39 rare missense variants normalized to WT (dotted line). The data are presented as percentages of WT ± S.D. (n = 2–4 independent transfections), with all points plotted. D, quantitation of expression levels of GPR39 rare SYN_ΔCB variants normalized to WT (dotted line). The details are the same as in C. Statistical analysis for both figures was two-tailed t test relative to WT. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.
      GPR39 signaling is complex, with multiple signaling outputs. However, in both endogenous and heterologous expression systems, GPR39 activates inositol trisphosphate (IP3)–dependent release of intracellular Ca2+ (
      • Asraf H.
      • Salomon S.
      • Nevo A.
      • Sekler I.
      • Mayer D.
      • Hershfinkel M.
      The ZnR/GPR39 interacts with the CaSR to enhance signaling in prostate and salivary epithelia.
      ,
      • Cohen L.
      • Asraf H.
      • Sekler I.
      • Hershfinkel M.
      Extracellular pH regulates zinc signaling via an Asp reside of the zinc-sensing receptor (ZnR/GPR39).
      ,
      • Verhulst P.J.
      • Lintermans A.
      • Janssen S.
      • Loeckx D.
      • Himmelreich U.
      • Buyse J.
      • Tack J.
      • Depoortere I.
      GPR39, a receptor of the ghrelin receptor family, plays a role in the regulation of glucose homeostasis in a mouse model of early onset diet-induced obesity.
      ). The first step in the signaling pathway is Gq-mediated activation of phospholipase C, which produces diacylglycerol and IP3. Here we exposed HEK293 cells transiently expressing WT or variant GPR39 to ZnCl2 and used a FRET assay to quantify accumulation of inositol monophosphate (IP1) in the presence of LiCl. Fig. 6A illustrates normalized maximal levels of IP1 for WT and representative missense variants having decreased (R386C) or increased (W314R) EC50 for activation by ZnCl2. Fig. 6C plots the EC50 values and 95% confidence intervals (CIs) for WT and all other missense variants, color-coded for their locations in GPR39 domains. Three variants, S78L, T122M, and R133G, had no measurable activity, likely accounted for by low levels of expression (recall Fig. 5, A and C), whereas M87I was well-expressed but inactive. One variant, R193C, had a dose-response curve that was significantly right-shifted, precluding determination of EC50 or Vmax, despite expression levels not significantly different from WT. Likewise, Fig. 6B illustrates the normalized dose-response relations for WT GPR39 and representative SYN_ΔCB variants that decrease (Val-354) or increase (Ile-29) the EC50 for ZnCl2 activation, and Fig. 6F plots the EC50 values and 95% CI of WT GPR39 and all SYN_ΔCB variants. In contrast to the missense variants, all SYN_ΔCB variants had ZnCl2-stimulated IP1 activity. Eight variants had increased EC50 values, and two variants had significantly reduced EC50 values relative to WT. Surface expression of all rare missense and SYN_ΔCB variants was assessed by BgTx-Alexa594 binding on nonpermeabilized cells (Fig. 6, E and H), and generally correlated with the extrapolated Vmax values for the IP1 assays (Fig. 6, D and G). The data were corrected for the very low levels of unspecific BgTx–Alexa 594 binding to untransfected HEK293 cells (6.03 ± 1.8% of WT (n = 6 independent experiments)). Total (Western blotting) and plasma membrane localization (BgTx binding) were comparable for many of the missense variants (H19R, A167S, S222Y, A293P, R302Q, P310A, R320W, R362C, and R386C), but total expression was greater than plasma membrane localization for five variants (M87I, C191W, R193C, A307V, and S336R), suggesting potential alterations in trafficking, whereas the remaining 12 variants had reduced total expression relative to plasma membrane localization, suggesting a potential increase in quality control-mediated degradation at the endoplasmic reticulum.
      Figure thumbnail gr6
      Figure 6Functional analysis of rare missense and synonymous GPR39 variants. A, IP1 responses to various concentration of ZnCl2 for all variants that expressed ≥50% of WT. Representative curves and fits for WT, W314R, and R386C are illustrated. Plotted are the means ± S.D. of three independent experiments, and lines show fits to Michaelis–Menten equation. The dotted line indicates 95% CI of WT. B, IP1 responses of representative SYN_ΔCB variants showing altered EC50 compared with WT, Ile-29, and Val-354. The details are as in A. C, EC50 values for dose-response relations of all functional rare missense variants (as in A), plotted as means ± 95% CI for n = 3 independent experiments. The dotted lines indicate the 95% CI range for WT. D, extrapolated Vmax of Michaelis–Menten fits to IP1 dose-response relations of rare missense variants, plotted as mean of three independent experiments plus 95% CI. The dotted lines indicate the 95% CI range for WT. E, bungarotoxin binding was used to assess surface localization of rare missense variants, plotted as mean plus all points, normalized to WT (dotted line), and corrected for unspecific binding in untransfected HEK293 cells. F, EC50 values of all rare SYN_ΔCB variants with measurable activity, analyzed and plotted as in C. G, extrapolated Vmax of fits to IP1 dose-response relations for rare SYN_ΔCB variants, analyzed, and plotted as in D. H, bungarotoxin binding to assess surface localization of rare SYN_ΔCB variants, analyzed, and plotted as in E. For all plots. *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001, by two-tailed t test relative to WT, assuming unequal variances.
      SYN_ΔCB plasma membrane localization varied more widely relative to WT, with a subset of variants having significantly elevated plasma membrane localization (Fig. 6H), which was not a reflection of increased expression assessed by Western blotting (Fig. 5D). Six variants had proportional total and plasma membrane localization (Ser-78, Ala-321, Ser-248, Leu-363, Pro-391, and Ala-394), whereas three had higher total expression versus plasma membrane localization (Ser-147, Gln-265, and Gly-4370), and the remainder had reduced total expression but elevated plasma membrane targeting. The most striking SYN_ΔCB variant, Gln-265, had an increased EC50 value relative to WT (Fig. 6F), increased Vmax, and surface localization relative to WT (Fig. 6, G and H), as well as increased total expression (Fig. 5D). Comparing total and plasma membrane expression of rare missense and SYN_ΔCB variants thus identified a subset of high-value variants that likely have significantly altered biosynthesis and/or trafficking, which could be confirmed by mechanistic studies.
      Overall, the results demonstrate a striking parallel in functional impact on expression and IP1 signaling for a subset of both rare missense and SYN_ΔCB variants of GPR39, rationalizing the disease associations identified for both classes of variants. The congruence of top disease associations between rare LOF, pathogenic missense, and SYN_ΔCB variants for a significant subset of GPCRs (recall Table 1) argues strongly for the utility of these analytical approaches for identification of functionally impactful rare variants and the imperative to include rare synonymous variants in genomic and functional analyses.

      Discussion

      The related GWAS and PheWAS approaches, which utilize common variants, have produced strong, validated targets for drug development (
      • Roth B.L.
      • Kroeze W.K.
      Integrated approaches for genome-wide interrogation of the druggable non-olfactory G protein–coupled receptor family.
      ). However, many genes do not have functionally important common variants, and additional approaches are required. The advantages of rare coding variant approaches are 2-fold. First, rare variants can have stronger effects on protein function, and second, assessment of the impact of rare coding variants on protein expression and function is straightforward. However, because of their low frequency, it is necessary to combine rare coding variants to increase the statistical power for detecting disease associations. In this study, we used pathogenicity triage as a basis for rare variant binning, based on the hypothesis that aggregating variants with similar directions of effect on protein function would amplify disease associations. We considered, in order of the strength of their potential effects, predicted LOF variants, predicted likely pathogenic or pathogenic missense variants, and synonymous variants with altered codon frequency bias. As a test of the importance of rare variant triage on disease associations, we compared agnostic binning of all rare variants in GPR39 with the disease associations identified for the three sorted variant classes. Of note, binning of all rare variants yielded no significant disease associations, whereas curated classes of those same rare variants had concordant top phenome-wide disease associations, obtained with nonoverlapping individuals. Likewise, SKAT analysis of the subset of rare missense variants predicted to be benign yielded no significant disease associations. The results argue strongly that binning all rare variants in a gene or genetic region dilutes the effects of impactful variants, and appropriate pathogenicity triage of rare variants is an integral first step in successful application of disease association analyses.
      Over the past 10 years, rare synonymous variants have been shown to impact protein expression and function with a frequency comparable with nonsynonymous variants, although the mechanism(s) by which they exert their biological function are diverse (
      • Hunt R.C.
      • Simhadri V.L.
      • Iandoli M.
      • Sauna Z.E.
      • Kimchi-Sarfaty C.
      Exposing synonymous mutations.
      • Bali V.
      • Bebok Z.
      Decoding mechanisms by which silent codon changes influence protein biogenesis and function.
      ,
      • Brar G.A.
      Beyond the triplet code: context cues transform translation.
      ,
      • Supek F.
      • Miñana B.
      • Valcárcel J.
      • Gabaldón T.
      • Lehner B.
      Synonymous mutations frequently act as driver mutations in human cancers.
      ,
      • Fernandez-Calero T.
      • Cabrera-Cabrera F.
      • Ehrlich R.
      • Marin M.
      Silent polymorphisms: can the tRNA population explain changes in protein properties?.
      ,
      • Jung H.
      • Lee D.
      • Lee J.
      • Park D.
      • Kim Y.J.
      • Park W.Y.
      • Hong D.
      • Park P.J.
      • Lee E.
      Intron retention is a widespread mechanism of tumor-suppressor inactivation.
      ,
      • Soussi T.
      • Taschner P.E.
      • Samuels Y.
      Synonymous somatic variants in human cancer are not infamous: a plea for full disclosure in databases and publications.
      ,
      • Kukreti R.
      • Tripathi S.
      • Bhatnagar P.
      • Gupta S.
      • Chauhan C.
      • Kubendran S.
      • Janardhan Reddy Y.C.
      • Jain S.
      • Brahmachari S.K.
      Association of DRD2 gene variant with schizophrenia.
      ,
      • Grymek K.
      • Lukasiewica S.
      • Faron-Góreckaa A.
      • Tworzydlo M.
      • Polit A.
      • Dziedzicka-Wasylewska M.
      Role of silent polymorphisms within the dopamine D1 receptor associated with schizophrenia on D1-D2 receptor heterodimerization.
      ,
      • Duan J.
      • Wainwright M.S.
      • Comeron J.M.
      • Saitou N.
      • Sanders A.R.
      • Gelernter J.
      • Gejman P.V.
      Synonymous mutations in the human dopamine receptor D2 (DRD2) affect mRNA stability and synthesis of the receptor.
      ,
      • Kimchi-Sarfaty C.
      • Oh J.M.
      • Kim I.-W.
      • Sauna Z.E.
      • Calcagno A.M.
      • Ambudkar S.V.
      • Gottesman M.M.
      A “silent” polymorphism in the MDR1 gene changes substrate specificity.
      • Shah K.
      • Cheng Y.
      • Hahn B.
      • Bridges R.
      • Bradbury N.A.
      • Mueller D.M.
      Synonymous codon usage affects the expression of wild type and F508del CFTR.
      ). Rare synonymous variants have not, however, been systematically utilized for disease association analysis (
      • Hunt R.C.
      • Simhadri V.L.
      • Iandoli M.
      • Sauna Z.E.
      • Kimchi-Sarfaty C.
      Exposing synonymous mutations.
      ,
      • Bali V.
      • Bebok Z.
      Decoding mechanisms by which silent codon changes influence protein biogenesis and function.
      ,
      • Soussi T.
      • Taschner P.E.
      • Samuels Y.
      Synonymous somatic variants in human cancer are not infamous: a plea for full disclosure in databases and publications.
      ,
      • Zhang T.
      • Wu Y.
      • Lan Z.
      • Shi Q.
      • Yang Y.
      • Guo J.
      Syntool: a novel region-based intolerance score to single nucleotide substitution for synonymous mutations predictions based on 12,136 individuals.
      ). Here we demonstrate that 39% of understudied or orphan GPCRs have impactful rare synonymous variants that yield disease associations comparable with LOF or missense variants in the same gene, which may provide impetus for exploratory drug development despite the lack of identified endogenous ligands. Our data also highlight the necessity for improved approaches for pathogenicity triage of rare synonymous variants. For GPR39, SKAT analysis on all rare synonymous variants for the most part yielded disease associations similar to those obtained with the SYN_ΔCB subset but have reduced significance, arguing for dilution of impactful variants. However, a few disease associations, including psoriatic arthropathy, showed more significant p values in the combined synonymous variant analysis (3.87E-11) versus SYN_ΔCB (1.03E-06), implying that simply isolating those synonymous variants having altered local codon usage bias missed some variants with significant functional impacts (see Table 2, Table 3). Using both LOF/missense and synonymous variants as independent data sets in the SKAT analyses increases confidence in the identified disease associations and provided the impetus for further functional studies.
      A disease association does not prove causation but implies an impact on protein expression and/or function that must be validated experimentally. As proof of concept, GPR39 was chosen because its agonist and dominant signaling pathway were known, and rare variant classes produced strongly congruent disease associations. We identified 168 rare variants in GPR39 in the 51,289 WES DiscovEHR cohort. Most individuals were heterozygous for the rare variants studied, and only 1 of 327 individuals had more than one rare variant. Of the 110 missense variants, the RMPath pipeline predicted that 28 were likely to be pathogenic (pLP + pP). Among these, rare variants previously shown to be important for GPR39 function were identified, including the zinc-binding site (
      • Storjohann L.
      • Holst B.
      • Schwartz T.W.
      Molecular mechanism of Zn2+ agonism in the extracellular domain of GPR39.
      ), a critical extracellular disulfide bond (
      • Storjohann L.
      • Holst B.
      • Schwartz T.W.
      A second disulfide bridge from the N-terminal domain to extracellular loop 2 dampens receptor activity in GPR39.
      ), the pH sensor (
      • Cohen L.
      • Asraf H.
      • Sekler I.
      • Hershfinkel M.
      Extracellular pH regulates zinc signaling via an Asp reside of the zinc-sensing receptor (ZnR/GPR39).
      ), and the NPXXY motif common to all family A GPCRs (
      • Yuan S.
      • Filipek S.
      • Palczewski K.
      • Vogel H.
      Activation of G-protein-coupled receptors correlates with the formation of a continuous internal water pathway.
      ). When expressed in HEK293 cells, these variants had the expected effects on EC50 values and maximal activities. Of the 47 synonymous variants, the codon usage table identified 21 variants having significantly altered codon usage frequency, and many had significant effects on expression, surface localization, and Zn2+-induced IP1 production. Surprisingly, the distributions of rare missense and SYN_ΔCB variants in the GPR39 structural domains were distinct. The majority of missense variants were localized to the transmembrane domain (77%), half in the extracellular loops (10 of 26), one an intracellular loop, and the remainder in the transmembrane helices (9 of 26). In contrast, ∼50% of SYN_ΔCB variants were in the C terminus; the remainder were distributed among transmembrane helices (3 of 20), intracellular (2 of 20), and extracellular (4 of 20) loops. The C-terminal variants in both classes had significant effects on expression and signaling, primarily reducing Vmax. In relation to many family A GPCRs, the C terminus of GPR39 is relatively large (108 amino acids) and has 22 putative phosphorylation sites. GPR39 signaling outputs are regulated by agonist-toggled binding of protein kinase inhibitor β (
      • Kovacs Z.
      • Schacht T.
      • Herrmann A.K.
      • Albrecht P.
      • Lefkimmiatis K.
      • Methner A.
      Protein kinase inhibitor β enhances the constitutive activity of G-protein-coupled zinc receptor GPR39.
      ). In addition, phosphorylation by Rho kinase mediates desensitization of GPR39 (
      • Shimizu Y.
      • Koyama R.
      • Kawamoto T.
      Rho kinase-dependent desensitization of GPR39: a unique mechanism of GPCR downregulation.
      ). Thus, C-terminal missense or synonymous variants may alter secondary structure and/or phosphorylation to impact signaling. The results focus attention on the C terminus as a critical domain mediating distinct aspects of GPR39 signaling. It should be noted that GPR39 signaling is complex, including Gs-mediated cAMP generation, Gq-mediated IP3 production, and serum response factor-serum response element (SRF-SRE)–mediated transcriptional regulation via G12/13 (
      • Shimizu Y.
      • Koyama R.
      • Kawamoto T.
      Rho kinase-dependent desensitization of GPR39: a unique mechanism of GPCR downregulation.
      ). The current study examined only Gq-mediated IP3 generation, leaving open the possible impact of missense and/or synonymous variants on signaling bias. Further, our results suggest that structural studies on a subset of rare missense and SYN_ΔCB variants having strong effects on signaling may reveal important structural determinants for optimal GPR39 function.
      Based on ClinVar and related databases, GPR39 has not been causally linked to disease. The present study focused on the general patient population, and all individuals in the present study were heterozygous for GPR39 rare variants. Nevertheless, both LOF/MISS and SYN_ΔCB variants yielded disease associations by SKAT, which suggests a role for GPR39 in nerve function, benign hyperplasia of the prostate, psoriatic arthropathy, diseases of hair and hair follicles, and benign essential hypertension. These associations are in agreement with cellular and animal studies suggesting that GPR39 regulates neuronal transport of potassium and chloride, attenuating postsynaptic excitability (
      • Sunuwar L.
      • Gilad D.
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ,
      • Chorin E.
      • Vinograd O.
      • Fleidervish I.
      • Gilad D.
      • Herrmann S.
      • Sekler I.
      • Aizenman E.
      • Hershfinkel M.
      Upregulation of KCC2 activity by zinc-mediated neurotransmission via the mZnR/GPR39 receptor.
      ,
      • Saadi R.A.
      • He K.
      • Hartnett K.A.
      • Kandler K.
      • Hershfinkel M.
      • Aizenman E.
      SNARE-dependent upregulation of potassium chloride co-transporter 2 activity after metabotropic zinc receptor activation in rat cortical neurons in vitro.
      • Gilad D.
      • Shorer S.
      • Ketzef M.
      • Friedman A.
      • Sekler I.
      • Aizenman E.
      • Hershfinkel M.
      Homeostatic regulation of KCC2 activity by the zinc receptor mZnR/GPR39 during seizures.
      ). Animal models demonstrate GPR39-mediated up-regulation of the neuronal K+/Cl cotransporter attenuates seizure activity (
      • Gilad D.
      • Shorer S.
      • Ketzef M.
      • Friedman A.
      • Sekler I.
      • Aizenman E.
      • Hershfinkel M.
      Homeostatic regulation of KCC2 activity by the zinc receptor mZnR/GPR39 during seizures.
      ), and knockout of GPR39 promotes depression (
      • Mlyniec K.
      • Budziszewska B.
      • Holst B.
      • Ostachowicz B.
      • Nowak G.
      GPR39 (zinc receptor) knockout mice exhibit depression-like behavior and CREB/BDNF down-regulation in the hippocampus.
      ). GPR39 is up-regulated by some anti-depressants (
      • Mlyniec K.
      • Nowak G.
      Up-regulation of the GPR39 Zn2+-sensing receptor and CREB/BDNF/TrkB pathway after chronic but not acute antidepressant treatment in the frontal cortex of zinc-deficient mice.
      ) in mouse models. One study in humans showed that GPR39 levels were significantly reduced in the hippocampus and cortex of suicide victims (
      • Mlyniec K.
      • Doboszewska U.
      • Szewczyk B.
      • Sowa-Kućma M.
      • Misztak P.
      • Piekoszewski W.
      • Trela F.
      • Ostachowicz B.
      • Nowak G.
      The involvement of the GPR39-Zn2+–sensing receptor in the pathophysiology of depression: studies in rodent models and suicide victims.
      ). The highest levels of Zn2+ in the body are in the prostate and seminal fluid (
      • Costello L.C.
      • Franklin R.B.
      A comprehensive review of the role of zinc in normal prostate function and metabolism; and its implications in prostate cancer.
      ), and GPR39 is differentially expressed in normal versus malignant prostate cells, where it regulates cell growth and proliferation (
      • Asraf H.
      • Salomon S.
      • Nevo A.
      • Sekler I.
      • Mayer D.
      • Hershfinkel M.
      The ZnR/GPR39 interacts with the CaSR to enhance signaling in prostate and salivary epithelia.
      ,
      • Sunuwar L.
      • Gilad D.
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ,
      • Hershfinkel M.
      The zinc sensing receptor, ZnR/GPR39, in health and disease.
      ), potentially contributing to the association with benign prostate hyperplasia observed in the present work. Finally, the observation that individuals heterozygous for rare variants in GPR39 have a significant association with diseases of hair and hair follicles should be considered in light of a recent study that found that GPR39 marks a class of stem cells in the sebaceous gland and contributes to wound healing in mice (
      • Zhao H.
      • Qiao J.
      • Zhang S.
      • Zhang H.
      • Lei H.
      • Lei X.
      • Wang X.
      • Deng Z.
      • Ning L.
      • Cao Y.
      • Guo Y.
      • Liu S.
      • Duan E.
      GPR39 marks specific cells within the sebaceous gland and contributes to skin wound healing.
      ). Overall, the identification of a subset of clinical phenotypes that echo cell culture and animal model studies supports disparate roles for GPR39 in human physiology and validates the SKAT approaches used herein to identify disease associations. Additional phenotypes, not previously associated with GPR39 including psoriatic arthropathy and macular puckering of the retina, deserve further attention.
      In conclusion, we have identified strong disease associations for GPR39 using either LOF/missense or rare synonymous variants that confirm and extend cellular and animal studies. Rare variants from both the missense and synonymous classes altered GPR39 expression and function, providing a rationale for the observed disease associations. The approach applied to a larger set of orphan or understudied GPCRs also yielded robust associations with phenome-wide significance, providing potential novel targets for further investigation. The methods described in this report are easily applied to any set of genes of potential clinical importance with available genomic and integrated longitudinal electronic health records. Finally, this study highlights the importance of rare synonymous variants in human physiology and argues for their routine inclusion in any comprehensive analysis of genomic variants as potential causes of disease.

      Experimental procedures

      The MyCode® Community Health Initiative recruits Geisinger patients from a broad range of inpatient and outpatient clinics to integrate genetic information with their clinical EHR data to foster discoveries in health and disease (
      • Carey D.J.
      • Fetterolf S.N.
      • Davis F.D.
      • Faucett W.A.
      • Kirchner H.L.
      • Mirshahi U.
      • Murray M.F.
      • Smelser D.T.
      • Gerhard G.S.
      • Ledbetter D.H.
      The Geisinger MyCode community health initiative: an electronic health record-linked biobank for precision medicine research.
      ). This study used only deidentified data and does not therefore involve “human subjects” as defined in 45CFR46.102(f). Although this study was not subject to oversight by the Geisinger Institutional Review Board, the specific uses of the genetic and associated clinical data described herein were approved by the Geisinger Institutional Review Board. We used WES from 51,289 participants (
      • Dewey F.E.
      • Murray M.F.
      • Overton J.D.
      • Habegger L.
      • Leader J.B.
      • Fetterolf S.N.
      • O'Dushlaine C.
      • Van Hout C.V.
      • Staples J.
      • Gonzaga-Jauregui C.
      • Metpally R.
      • Pendergrass S.A.
      • Giovanni M.A.
      • Kirchner H.L.
      • Balasubramanian S.
      • et al.
      Distribution and clinical impact of functional variants in 50,726 whole-exome sequences from the DiscovEHR study.
      ), which includes patients enrolled through primary care and specialty outpatient clinics. Participants were 59% female, with a median age of 61 years, and predominantly Caucasian (98%) and non-Hispanic/Latino (95%). Patient demographics are in Table 4. Whole exome sequencing was performed using 75-bp paired-end sequencing on an Illumina v4 HiSeq 2500 to a coverage depth sufficient to provide greater than 20× haploid read depth of over 85% of targeted bases for at least 96% of samples. Sample preparation, sequencing, sequence alignment, variant identification, genotype assignment, and quality control steps were carried out as previously described (
      • Carey D.J.
      • Fetterolf S.N.
      • Davis F.D.
      • Faucett W.A.
      • Kirchner H.L.
      • Mirshahi U.
      • Murray M.F.
      • Smelser D.T.
      • Gerhard G.S.
      • Ledbetter D.H.
      The Geisinger MyCode community health initiative: an electronic health record-linked biobank for precision medicine research.
      ).
      Table 4Demographics of the DiscovEHR cohort used in the present study
      Basic demographicsDiscoveEHR sequenced patients
      Total patients, N51,289
      Female, N (%)30,290 (59%)
      Median age, years61 (range 48–72)
      Median body mass index, kg/m230 (range 26–36)
      Median years of electronic health record data14 (range 11–17)
      Median medication orders/patient129 (range 37–221)
      Rare variants were separated into two groups: 1) nonsense and frameshift variants and 2) missense variants. Nonsense and frameshift variants were scored by SnpEff fields LOF[*]·PERC (percent) ≥ 0.5 or (nonsense mediated decay NMD)[*]. Number of transcripts (NUMTR) values > 3 are considered potential LOF variants. To score missense variants, we summed numerical values derived from pathogenicity prediction tools to generate a cumulative score as follows (RMPath, in silico radical mutations pathogenicity prediction; maximum score = 24): 1) SIFT scoring: D (damaging) = 2; all others = 0; 2) PolyPhen2 scoring: (i) HDIV: probably damaging (D) = 2; possibly damaging (P) = 1; all other results = 0. (ii) HVAR: probably damaging (D) = 2; possibly damaging (P) = 1; all other results = 0; 3) MutationTaster scoring: Disease causing (D) = 2; all other results = 0; 4) MutationAssessor scoring: high (H) = 2; medium (M) = 1; all other results = 0; 5) FATHMM scoring: D (damaging) = 2; all other results = 0; 6) LRT: D (deleterious) → score 2; All others → score 0; 7) PROVEAN scoring: D (damaging) = 2; all other results = 0; 8) GRANTHAM scoring: ≥140 = 4; ≥120 = 3; ≥90 = 2; ≥70 = 1; all other results = 0; and 9) Genomic Evolutionary Rate Profiling GERP++ rejected substitutions (RS) scoring: ≥6 = 4; ≥5 = 3; ≥3.5 = 2; ≥2 = 1; all other results = 0. Variants with RMPath scores ≤8 were pB, those with scores in the range from 6 to 11 were predicted to be Likely Benign (pLB); those with scores between 11 and 17 were predicted to be likely pathogenic (pLP); and those with scores ≥17 were predicted to be pathogenic (pP). For clinical purposes, all variants are considered variants of unknown significance until family cosegregation and/or functional analyses definitively establish a variant as pathogenic. Our goal in this study was to identify variants likely to impact GPR39 expression or function.
      For all rare synonymous codons, we determined the ratio of the genome-wide codon fraction of the GPR39 reference codon over the synonymous codon variant at that position using human-specific codon usage tables (
      • Nakamura Y.
      • Gojobori T.
      • Ikemura T.
      Codon usage tabulated from international DNA sequence databases: status for the year 2000.
      ,
      • Athey J.
      • Alexaki A.
      • Osipova E.
      • Rostovtsev A.
      • Santana-Quintero L.V.
      • Katneni U.
      • Simonyan V.
      • Kimchi-Sarfaty C.
      A new and updated resource for codon usage tables.
      ), e.g. for S147S (tcG/tcA), 0.15/0.06 = 2.5, i.e. a 2.5-fold change. All synonymous variants with at least a 2-fold change in site-specific codon usage were used in the SKAT analyses.
      Deidentified EHR data were obtained from an approved data broker. All unique ICD9 codes with ≥200 patients (regardless of genotype) with at least three independent incidences of the ICD9 code were extracted from the EHR (
      • Verma A.
      • Bradford Y.
      • Dudek S.
      • Lucas A.M.
      • Verma S.S.
      • Pendergrass S.A.
      • Ritchie M.D.
      A simulation study investigating power estimates in phenome-wide association studies.
      ). Individuals having one or two calls were excluded, and those having no calls of a particular ICD9 code were considered as controls. All non-Europeans were excluded from the analysis, as was one sample from pairs of closely related individuals up to first cousins. All models were adjusted for sex, age, age2, and the first four principal components. The PheWAS R package was used for association analyses, which were performed for three GPR39 common variants (2:133174764, pA50V; 2:133174999, pT128T; and 2:133402607, noncoding), and plotting of results was done using GraphPad Prism (version 6). We used the sequence kernel association test with default weights (SKAT R package) to compare the burden of rare variants in cases and controls. Cases and controls were defined as for PheWAS analysis (
      • Verma A.
      • Bradford Y.
      • Dudek S.
      • Lucas A.M.
      • Verma S.S.
      • Pendergrass S.A.
      • Ritchie M.D.
      A simulation study investigating power estimates in phenome-wide association studies.
      ). Age, age2, sex, body mass index, and the first four principal components were used as covariates. All non-European samples were excluded. We used principal components calculated from high quality common SNPs to estimate the ancestry of our samples. Using a kernel density estimator, the likelihood that an individual belonged to each HapMap class (African, admixed American, East Asian, European, or south Asian) was estimated; samples that fell into an unknown group were excluded (
      • Staples J.
      • Maxwell E.K.
      • Gosalia N.
      • Gonzaga-Jauregui C.
      • Snyder C.
      • Hawes A.
      • Penn J.
      • Ulloa R.
      • Bai X.
      • Lopez A.E.
      • Van Hout C.V.
      • O'Dushlaine C.
      • Teslovich T.M.
      • McCarthy S.E.
      • Balasubramanian S.
      • et al.
      Profiling and leveraging relatedness in a Precision Medicine Cohort of 92,455 exomes.
      ). Statistical significance was determined using the Bonferroni correction, α/m, where α is 0.05, and m is the number of variant sets and disease codes, as indicated for each analysis (
      • Verma A.
      • Lucas A.
      • Verma S.S.
      • Zhang Y.
      • Josyula N.
      • Khan A.
      • Hartzel D.N.
      • Lavage D.R.
      • Leader J.
      • Ritchie M.D.
      • Pendergrass S.A.
      PheWAS and beyond: the landscape of associations with medical diagnoses and clinical measures across 38,662 individuals from Geisinger.
      ).
      All missense and SYN_ΔCB variants were generated by gene synthesis (Thermo Fisher Scientific, Gene Art) based on the NP_001499.1 reference sequence (GeneID 2863, 453 amino acids) and included tandem tags at the N terminus to facilitate expression and functional studies. After the initiation methionine codon, we added a 3× FLAG epitope plus a minimal bungarotoxin-binding site as shown: GACTACAAAGACCATGACGGTGATTATAAAGATCATGATATCGATTACAAGGATGACGATGACAAGGGTTGGAGATACTACGAGAGCTCCCTGGAGCCCTACCCTGACGGT.
      WT and variant constructs were transiently transfected into HEK293 cells. The cells were lysed after 48 h, and 25 μg of lysates were separated on 4–15% SDS-PAGE gels (Bio-Rad), blotted to nitrocellulose, blocked with 5% milk in TBS-T, and exposed for 1 h at room temperature to anti-FLAG mAb conjugated to HRP (Sigma). The blots were exposed to SuperSignal West Pico Chemiluminescence Substrate (Thermo Fisher), and the images were recorded on a FUJIFILM LAS-4000mini luminescence analyzer and processed with Image-Gauge version 3.0. GPR39 expression was quantified using Multi-Gauge software, normalized to WT expression, and corrected for HEK293 background on the same blot. Each variant was assayed in two to four independent transfections. To accurately quantify the range of expression observed for SYN_ΔCB variants, we segregated high- and low-expressing variants plus WT onto separate gels, allowing differential exposure times for chemiluminescence acquisition.
      WT and variant constructs were transiently transfected into HEK293 cells and cultured for 24 h, and equivalent numbers of cells were plated in poly-d-lysine–coated 96-well plates for a further 24 h. Replicates were treated with various concentrations of ZnCl2 for 10 min at 37 °C. Levels of IP1 were determined as recommended by manufacturer (IP-One ELISA kit, Cisbio USA, Inc.). The plates were read on a POLARStar Omega (BMG Labtech) plate reader. We made no assumptions regarding signaling mechanism and its potential cooperativity and estimated Vmax and EC50 from Zn2+ versus IP1 data combined from three independent experiments using GraphPad PRISM (V6) and a general dose-response equation with the potential for cooperativity (Y = Y0 + [(VmaxY0)/(1 + 10 (logEC50 − [Zn2+] n)]), where Y0 is baseline IP1 in the absence of Zn2+, Vmax is the maximal response, EC50 is the Zn2+ concentration at 50% of Vmax, and n is the Hill coefficient (slope). Parallel plates from three independent transfections were exposed to Alexa 594–conjugated bungarotoxin (3 μg/ml, 5 min at 37 °C, 5% CO2) to quantify surface expression of GPR39 variants. After removal of bungarotoxin and washes with PBS, fluorescence was read on the POLARStar Omega (BMG Labtech) plate reader. WT, standards, and solution blanks were included in each assay, and unspecific binding of BgTx–Alexa 594 to HEK293 cells was determined relative to cells transfected with WT GPR39 in the same experiment.

      Author contributions

      R. D., R. P. R. M., K. J., S. K., and G. E. B. investigation; R. D., R. P. R. M., K. J., S. K., D. T. S., R. G. C., and G. E. B. visualization; R. D. and R. G. C. methodology; R. D., R. P. R. M., D. T. S., M. H., D. J. C., J. D. R., and G. E. B. writing-review and editing; R. P. R. M., R.G.C., and G. E. B. data curation; R. P. R. M. software; R. P. R. M., K. J., D. T. S., and G. E. B. formal analysis; R. P. R. M. and G. E. B. supervision; D. T. S. validation; R. G. C. resources; G. E. B. conceptualization; G. E. B. funding acquisition; G. E. B. writing-original draft; G. E. B. project administration.

      Acknowledgments

      We are grateful to Geisinger and the Regeneron Genetics Center (G. Abecasis, A. Baras, M. Cantor, G. Coppola, A. Economides, J.D. Overton, J. G. Reid, A. Shuldiner, C. Beechert, C. Forsythe, E. D. Fuller, Z. Gu, M. Lattari, A. Lopez, T. D. Schleicher, M. S. Padilla, K. Toledo, L. Widom, S. E. Wolf, M. Pradhan, K. Manoochehri, R. H. Ulloa, X. Bai, S. Balasubramanian, L. Barnard, A. Blumenfeld, Y. Chai, G. Eom, L. Habegger, Y. Hahn, A. Hawes, S. Khalid, E. K. Maxwell, J. Penn, J. C. Staples, A. Yadav, P. M. Guzzardo, M. B. Jones, and L. J. Mitnaul) for access to the DiscovEHR cohort and to the patients in the MyCode Biobank.

      References

        • Hauser A.S.
        • Attwood M.M.
        • Rask-Andersen M.
        • Schiöth H.B.
        • Gloriam D.E.
        Trends in GPCR drug discovery: new agents, targets and indications.
        Nat. Rev. Drug Discov. 2017; 16 (10.1038/nrd.2017.178, 29075003): 829-842
        • Roth B.L.
        • Kroeze W.K.
        Integrated approaches for genome-wide interrogation of the druggable non-olfactory G protein–coupled receptor family.
        J. Biol. Chem. 2017; 290 (10.1074/jbc.R115.654764, 26100629): 19471-19477
        • Kovacs P.
        • Schönberg T.
        The relevance of genomic signatures at adhesion GPCR loci in humans.
        Handb. Exp. Pharmacol. 2016; 234 (10.1007/978-3-319-41523-9_9, 27832489): 179-217
        • Wu M.C.
        • Lee S.
        • Cai T.
        • Li Y.
        • Boehnke M.
        • Lin X.
        Rare-variant association testing for sequencing data with the sequence kernel association test.
        Am. J. Hum. Genet. 2011; 89 (10.1016/j.ajhg.2011.05.029, 21737059): 82-93
        • Lee S.
        • Emond M.J.
        • Bamshad M.J.
        • Barnes K.C.
        • Rieder M.J.
        • Nickerson D.A.
        • NHLBI GO Exome Sequencing Project–ESP Lung Project Team
        • Christiani D.C.
        • Wurfel M.M.
        • Lin X.
        Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies.
        Am. J. Hum. Genet. 2012; 91 (10.1016/j.ajhg.2012.06.007, 22863193): 224-237
        • Richardson T.G.
        • Shihab H.A.
        • Rivas M.A.
        • McCarthy M.I.
        • Campbell C.
        • Timpson N.J.
        • Gaunt T.R.
        A protein domain and family based approach to rare variant association analysis.
        PLoS One. 2016; 11 (10.1371/journal.pone.0153803, 27128313): e0153803
        • Friedrichs S.
        • Malzahn D.
        • Pugh E.W.
        • Almeida M.
        • Liu X.Q.
        • Bailey J.N.
        Filtering genetic variants and placing informative priors based on putative biological function.
        BCM Genet. 2016; 17 (10.1186/s12863-015-0313-x, 26866982): 8
        • Sriram K.
        • Insel P.A.
        G protein–coupled receptors as targets for approved drugs: how many targets and how many drugs?.
        Mol. Pharmacol. 2018; 93 (10.1124/mol.117.111062, 29298813): 251-258
        • Dewey F.E.
        • Murray M.F.
        • Overton J.D.
        • Habegger L.
        • Leader J.B.
        • Fetterolf S.N.
        • O'Dushlaine C.
        • Van Hout C.V.
        • Staples J.
        • Gonzaga-Jauregui C.
        • Metpally R.
        • Pendergrass S.A.
        • Giovanni M.A.
        • Kirchner H.L.
        • Balasubramanian S.
        • et al.
        Distribution and clinical impact of functional variants in 50,726 whole-exome sequences from the DiscovEHR study.
        Science. 2016; 354 (10.1126/science.aaf6814, 28008009): aaf6814
        • Asraf H.
        • Salomon S.
        • Nevo A.
        • Sekler I.
        • Mayer D.
        • Hershfinkel M.
        The ZnR/GPR39 interacts with the CaSR to enhance signaling in prostate and salivary epithelia.
        J. Cell. Physiol. 2014; 229 (10.1002/jcp.24514, 24264723): 868-877
        • Cohen L.
        • Asraf H.
        • Sekler I.
        • Hershfinkel M.
        Extracellular pH regulates zinc signaling via an Asp reside of the zinc-sensing receptor (ZnR/GPR39).
        J. Biol. Chem. 2012; 287 (10.1074/jbc.M112.372441, 22879599): 33339-33350
        • Sunuwar L.
        • Gilad D.
        • Hershfinkel M.
        The zinc sensing receptor, ZnR/GPR39, in health and disease.
        Front. Biosci. (Landmark Ed). 2017; 22 (10.2741/4554, 28199213): 1469-1492
        • Hershfinkel M.
        The zinc sensing receptor, ZnR/GPR39, in health and disease.
        Int. J. Mol. Sci. 2018; 19 (10.3390/ijms19020437, 29389900): E439
        • Basile A.O.
        • Byrska-Bishop M.
        • Wallace J.
        • Frase A.T.
        • Ritchie M.D.
        Novel features and enhancements in BioBin, a tool for biologically inspired binning and association analysis of rare variants.
        Bioinformatics. 2018; 34 (10.1093/bioinformatics/btx559, 28968757): 527-529
        • Giddens M.M.
        • Wong J.C.
        • Schroeder J.P.
        • Farrow E.G.
        • Smith B.M.
        • Owino S.
        • Soden S.E.
        • Meyer R.C.
        • Saunders C.
        • LePichon J.B.
        • Weinshenker D.
        • Escayg A.
        • Hall R.A.
        GPR37L1 modulates seizure susceptibility: evidence from mouse studies and analyses of a human GPR37L1 variant.
        Neurobiol. Dis. 2017; 106 (10.1016/j.nbd.2017.07.006, 28688853): 181-190
        • Hoeck J.D.
        • Biehs B.
        • Kurtova A.V.
        • Kljavin N.M.
        • de Sousa E.
        • Melo F.
        • Alicke B.
        • Koeppen H.
        • Modrusan Z.
        • Piskol R.
        • de Sauvage F.J.
        Stem cell plasticity enables hair regeneration following LGR5+ cell loss.
        Nat. Cell Biol. 2017; 19 (10.1038/ncb3535, 28553937): 666-676
        • Michel L.
        • Reygagne P.
        • Benech P.
        • Jean-Louis F.
        • Scalvino S.
        • Ly Ka So S.
        • Hamidou Z.
        • Bianovici S.
        • Pouch J.
        • Ducos B.
        • Bonnet M.
        • Bensussan A.
        • Patatian A.
        • Lati E.
        • Wdzieczak-Bakala J.
        • et al.
        Study of gene expression alteration in male androgenic alopecia: evidence of predominant molecular signaling pathways.
        Brit. J. Dermatol. 2017; 177 (10.1111/bjd.15577): 1322-1336
        • Guillemot-Legris O.
        • Mutemberezi V.
        • Cani P.D.
        • Muccioli G.G.
        Obesity is associated with changes in oxysterol metabolism and levels in mice liver, hypothalamus, adipose tissue and plasma.
        Sci. Rep. 2016; 6 (10.1038/srep19694, 26795945): 19694
        • Li Y.
        • Dong J.
        • Xuan Y.-H.
        • Liu S.-S.
        • Luo J.-Y.
        • Xiao Y.-J.
        • Tian Z.P.
        • Sun Z.J.
        Identification of microRNA expression in a rat model of post-infarction heart failure.
        Int. J. Clin. Exp. Pathol. 2017; 10: 4729-4738
        • Storjohann L.
        • Holst B.
        • Schwartz T.W.
        A second disulfide bridge from the N-terminal domain to extracellular loop 2 dampens receptor activity in GPR39.
        Biochemistry. 2008; 47 (10.1021/bi8005016, 18693759): 9198-9207
        • Storjohann L.
        • Holst B.
        • Schwartz T.W.
        Molecular mechanism of Zn2+ agonism in the extracellular domain of GPR39.
        FEBS Lett. 2008; 582 (10.1016/j.febslet.2008.06.030, 18588883): 2583-2588
        • Yuan S.
        • Filipek S.
        • Palczewski K.
        • Vogel H.
        Activation of G-protein-coupled receptors correlates with the formation of a continuous internal water pathway.
        Nat. Commun. 2014; 5 (10.1038/ncomms5733, 25203160): 4733
        • Carey D.J.
        • Fetterolf S.N.
        • Davis F.D.
        • Faucett W.A.
        • Kirchner H.L.
        • Mirshahi U.
        • Murray M.F.
        • Smelser D.T.
        • Gerhard G.S.
        • Ledbetter D.H.
        The Geisinger MyCode community health initiative: an electronic health record-linked biobank for precision medicine research.
        Genet. Med. 2016; 18 (10.1038/gim.2015.187, 26866580): 906-913
        • Verhulst P.J.
        • Lintermans A.
        • Janssen S.
        • Loeckx D.
        • Himmelreich U.
        • Buyse J.
        • Tack J.
        • Depoortere I.
        GPR39, a receptor of the ghrelin receptor family, plays a role in the regulation of glucose homeostasis in a mouse model of early onset diet-induced obesity.
        J. Neuroendocrinol. 2011; 23 (10.1111/j.1365-2826.2011.02132.x, 21470317): 490-500
        • Hunt R.C.
        • Simhadri V.L.
        • Iandoli M.
        • Sauna Z.E.
        • Kimchi-Sarfaty C.
        Exposing synonymous mutations.
        Trends Genet. 2014; 30 (10.1016/j.tig.2014.04.006, 24954581): 308-321
        • Bali V.
        • Bebok Z.
        Decoding mechanisms by which silent codon changes influence protein biogenesis and function.
        Int. J. Biochem. Cell Biol. 2015; 64 (10.1016/j.biocel.2015.03.011, 25817479): 58-74
        • Brar G.A.
        Beyond the triplet code: context cues transform translation.
        Cell. 2016; 167 (10.1016/j.cell.2016.09.022, 27984720): 1681-1692
        • Supek F.
        • Miñana B.
        • Valcárcel J.
        • Gabaldón T.
        • Lehner B.
        Synonymous mutations frequently act as driver mutations in human cancers.
        Cell. 2014; 156 (10.1016/j.cell.2014.01.051, 24630730): 1324-1335
        • Fernandez-Calero T.
        • Cabrera-Cabrera F.
        • Ehrlich R.
        • Marin M.
        Silent polymorphisms: can the tRNA population explain changes in protein properties?.
        Life (Basel). 2016; 6 (10.3390/life6010009, 26901226): E9
        • Jung H.
        • Lee D.
        • Lee J.
        • Park D.
        • Kim Y.J.
        • Park W.Y.
        • Hong D.
        • Park P.J.
        • Lee E.
        Intron retention is a widespread mechanism of tumor-suppressor inactivation.
        Nat. Genet. 2015; 47 (10.1038/ng.3414, 26437032): 1242-1248
        • Soussi T.
        • Taschner P.E.
        • Samuels Y.
        Synonymous somatic variants in human cancer are not infamous: a plea for full disclosure in databases and publications.
        Hum. Mutat. 2017; 38 (10.1002/humu.23163, 28026089): 339-342
        • Kukreti R.
        • Tripathi S.
        • Bhatnagar P.
        • Gupta S.
        • Chauhan C.
        • Kubendran S.
        • Janardhan Reddy Y.C.
        • Jain S.
        • Brahmachari S.K.
        Association of DRD2 gene variant with schizophrenia.
        Neurosci. Lett. 2006; 392 (10.1016/j.neulet.2005.08.059, 16183199): 68-71
        • Grymek K.
        • Lukasiewica S.
        • Faron-Góreckaa A.
        • Tworzydlo M.
        • Polit A.
        • Dziedzicka-Wasylewska M.
        Role of silent polymorphisms within the dopamine D1 receptor associated with schizophrenia on D1-D2 receptor heterodimerization.
        Pharmacol. Rep. 2009; 61 (10.1016/S1734-1140(09)70164-1, 20081237): 1024-1033
        • Duan J.
        • Wainwright M.S.
        • Comeron J.M.
        • Saitou N.
        • Sanders A.R.
        • Gelernter J.
        • Gejman P.V.
        Synonymous mutations in the human dopamine receptor D2 (DRD2) affect mRNA stability and synthesis of the receptor.
        Hum. Mol. Genet. 2003; 12 (10.1093/hmg/ddg055, 12554675): 205-216
        • Kimchi-Sarfaty C.
        • Oh J.M.
        • Kim I.-W.
        • Sauna Z.E.
        • Calcagno A.M.
        • Ambudkar S.V.
        • Gottesman M.M.
        A “silent” polymorphism in the MDR1 gene changes substrate specificity.
        Science. 2007; 315 (10.1126/science.1135308, 17185560): 525-528
        • Shah K.
        • Cheng Y.
        • Hahn B.
        • Bridges R.
        • Bradbury N.A.
        • Mueller D.M.
        Synonymous codon usage affects the expression of wild type and F508del CFTR.
        J. Mol. Biol. 2015; 427 (10.1016/j.jmb.2015.02.003, 25676312): 1464-1479
        • Zhang T.
        • Wu Y.
        • Lan Z.
        • Shi Q.
        • Yang Y.
        • Guo J.
        Syntool: a novel region-based intolerance score to single nucleotide substitution for synonymous mutations predictions based on 12,136 individuals.
        BioMed Res. Int. 2017; 2017 (28812016): 5096208
        • Kovacs Z.
        • Schacht T.
        • Herrmann A.K.
        • Albrecht P.
        • Lefkimmiatis K.
        • Methner A.
        Protein kinase inhibitor β enhances the constitutive activity of G-protein-coupled zinc receptor GPR39.
        Biochem. J. 2014; 462 (10.1042/BJ20131198, 24869658): 125-132
        • Shimizu Y.
        • Koyama R.
        • Kawamoto T.
        Rho kinase-dependent desensitization of GPR39: a unique mechanism of GPCR downregulation.
        Biochem. Pharmacol. 2017; 140 (10.1016/j.bcp.2017.06.115, 28619258): 105-114
        • Chorin E.
        • Vinograd O.
        • Fleidervish I.
        • Gilad D.
        • Herrmann S.
        • Sekler I.
        • Aizenman E.
        • Hershfinkel M.
        Upregulation of KCC2 activity by zinc-mediated neurotransmission via the mZnR/GPR39 receptor.
        J. Neurosci. 2011; 31 (10.1523/JNEUROSCI.2205-11.2011, 21900570): 12916-12926
        • Saadi R.A.
        • He K.
        • Hartnett K.A.
        • Kandler K.
        • Hershfinkel M.
        • Aizenman E.
        SNARE-dependent upregulation of potassium chloride co-transporter 2 activity after metabotropic zinc receptor activation in rat cortical neurons in vitro.
        Neuroscience. 2012; 210 (10.1016/j.neuroscience.2012.03.001, 22441041): 38-46
        • Gilad D.
        • Shorer S.
        • Ketzef M.
        • Friedman A.
        • Sekler I.
        • Aizenman E.
        • Hershfinkel M.
        Homeostatic regulation of KCC2 activity by the zinc receptor mZnR/GPR39 during seizures.
        Neurobiol. Dis. 2015; 81 (10.1016/j.nbd.2014.12.020, 25562657): 4-13
        • Mlyniec K.
        • Budziszewska B.
        • Holst B.
        • Ostachowicz B.
        • Nowak G.
        GPR39 (zinc receptor) knockout mice exhibit depression-like behavior and CREB/BDNF down-regulation in the hippocampus.
        Int. J. Neuropsychopharmacol. 2014; 18 (25609596): pyu002
        • Mlyniec K.
        • Nowak G.
        Up-regulation of the GPR39 Zn2+-sensing receptor and CREB/BDNF/TrkB pathway after chronic but not acute antidepressant treatment in the frontal cortex of zinc-deficient mice.
        Pharmacol. Rep. 2015; 67 (10.1016/j.pharep.2015.04.003, 26481532): 1135-1140
        • Mlyniec K.
        • Doboszewska U.
        • Szewczyk B.
        • Sowa-Kućma M.
        • Misztak P.
        • Piekoszewski W.
        • Trela F.
        • Ostachowicz B.
        • Nowak G.
        The involvement of the GPR39-Zn2+–sensing receptor in the pathophysiology of depression: studies in rodent models and suicide victims.
        Neuropharmacology. 2014; 79 (10.1016/j.neuropharm.2013.12.001, 24333148): 290-297
        • Costello L.C.
        • Franklin R.B.
        A comprehensive review of the role of zinc in normal prostate function and metabolism; and its implications in prostate cancer.
        Arch. Biochem. Biophys. 2016; 611 (10.1016/j.abb.2016.04.014, 27132038): 100-112
        • Zhao H.
        • Qiao J.
        • Zhang S.
        • Zhang H.
        • Lei H.
        • Lei X.
        • Wang X.
        • Deng Z.
        • Ning L.
        • Cao Y.
        • Guo Y.
        • Liu S.
        • Duan E.
        GPR39 marks specific cells within the sebaceous gland and contributes to skin wound healing.
        Sci. Rep. 2015; 5 (10.1038/srep07913, 25604641): 7913
        • Nakamura Y.
        • Gojobori T.
        • Ikemura T.
        Codon usage tabulated from international DNA sequence databases: status for the year 2000.
        Nucleic Acids Res. 2000; 28 (10.1093/nar/28.1.292, 10592250): 292
        • Athey J.
        • Alexaki A.
        • Osipova E.
        • Rostovtsev A.
        • Santana-Quintero L.V.
        • Katneni U.
        • Simonyan V.
        • Kimchi-Sarfaty C.
        A new and updated resource for codon usage tables.
        BMC Bioinformatics. 2017; 18 (10.1186/s12859-017-1793-7, 28865429): 391
        • Verma A.
        • Bradford Y.
        • Dudek S.
        • Lucas A.M.
        • Verma S.S.
        • Pendergrass S.A.
        • Ritchie M.D.
        A simulation study investigating power estimates in phenome-wide association studies.
        BMC Bioinformatics. 2018; 19 (10.1186/s12859-018-2135-0, 29618318): 120
        • Staples J.
        • Maxwell E.K.
        • Gosalia N.
        • Gonzaga-Jauregui C.
        • Snyder C.
        • Hawes A.
        • Penn J.
        • Ulloa R.
        • Bai X.
        • Lopez A.E.
        • Van Hout C.V.
        • O'Dushlaine C.
        • Teslovich T.M.
        • McCarthy S.E.
        • Balasubramanian S.
        • et al.
        Profiling and leveraging relatedness in a Precision Medicine Cohort of 92,455 exomes.
        Am. J. Hum. Genet. 2018; 102 (10.1016/j.ajhg.2018.03.012, 29727688): 874-889
        • Verma A.
        • Lucas A.
        • Verma S.S.
        • Zhang Y.
        • Josyula N.
        • Khan A.
        • Hartzel D.N.
        • Lavage D.R.
        • Leader J.
        • Ritchie M.D.
        • Pendergrass S.A.
        PheWAS and beyond: the landscape of associations with medical diagnoses and clinical measures across 38,662 individuals from Geisinger.
        Am. J. Hum. Genet. 2018; 102 (10.1016/j.ajhg.2018.02.017, 29606303): 592-608