Sex differences in the genetic architecture of obsessive-compulsive disorder.

The datasets ( Supplementary Figure 1 ) used in this study comprise the OCD Psychiatric Genomics Consortium sample and are fully described in primary publications ( Stewart et al. 2013 ; Mattheisen et al. 2015 ; International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS) 2017 ). All participants over 18 and the parents of participants under 18 gave written informed consent and this work was approved by the relevant institutional review boards at all participating sites. Participants of European ancestry were selected for this study and include cases and controls from Dutch, South African, European, and Ashkenazi Jewish ancestries. Additionally, trio samples were included in the meta-analysis and consisted of proband cases and pseudo-controls. The pseudo-controls were derived from the non-transmitted parental chromosomes. Genotype level data from all studies were pre-phased with SHAPEIT2 ( Delaneau, Zagury, and Marchini 2013 ), and imputed to the 1000 Genomes Project reference panel (Phase I integrated variant set release; NCBI build 37 (hg19)) using IMPUTE2 ( Howie, Marchini, and Stephens 2011 ), using the Ricopili pipeline ( Schizophrenia Working Group of the Psychiatric Genomics Consortium 2014 ). Prior to imputation, SNPs with call rate<0.98, minor allele frequency (MAF)0.02, Hardy–Weinberg equilibrium (HWE) p-values <1e-6 for controls and <1e-10 for cases were removed using PLINK ( Purcell et al. 2007 ). After imputation, any SNPs with IMPUTE2 info score <0.6 and certainty <0.8 were removed. After splitting the datasets by sex, SNPs with MAF <0.05 were removed from each sex, because the number of subjects carrying the minor allele was too small and could give rise to false positive association results. At the individual level, samples were removed if the genotyping call rate was 0.20, or there was an inconsistency between genetic sex and reported sex. Furthermore, pairwise identity by descent (IBD) analysis was used to identify cryptic relatedness between individuals, and one individual was removed at random from any pair related at the approximate level of first cousins (pi-hat>0.2). Principal component analyses were performed using EIGENSOFT ( Price et al. 2006 ) separately for each sub-population, ( Supplementary Methods ; Supplementary Figure 2 and 3 ) and PC plots were inspected to ensure that for every cluster of cases there was a cluster of controls in the same PC space. PCs are visualized separately by sex in order to ensure no subpopulation structure within each sex. After quality control, the total sample comprised 4,038 males and 5,832 females. The numbers of post-QC SNPs and individuals are listed in Supplementary Table 1 . X chromosome genotypes were processed separately from autosomal genotypes as additional care is required for pre-phasing, imputation, and post-imputation QC. At the genotype level, the pre-imputation QC steps for the X chromosome SNPs were the same as for the autosomes. An additional flag of -chrX was added when running SHAPEIT2 and IMPUTE2 software. Post-imputation, we employed the XWAS QC pipeline to remove variants in the pseudoautosomal regions (PARs), variants that were not in Hardy-Weinberg equilibrium in females, or variants with significantly different MAF (p<0.05/#SNPs) and differential missingness (P0.98 on the X chromosome. Furthermore, because we could not use the same case/pseudo-control design for the trio data (i.e. due to lack of a non-transmitted X chromosome from the fathers of affected females), we included only the affected individuals from the trio data, ancestry-matched them to controls from the case/control dataset, and analyzed them with the case/control data. We performed PCA using EIGENSOFT and removed any trio cases without matched controls. For each individual dataset, we performed sex-stratified and combined GWAS on imputed dosage files for autosomes and the X chromosome. In the sex-stratified analysis, dosages (i.e. number of chromosome copies) for the X chromosome in cases are equivalent to controls within each sex. However, in the sex-combined analysis, differences in dosage compensation between males (with one ChrX) and females (with two ChrX) should be considered. Thus, we verified that performing association analysis on ChrX dosage files produced consistent results compared with analysis of best-guess data in which ChrX was coded as 0/2 for males instead of 0/1 (--xchr-model 2 in PLINK). In all association analyses, principal components correlated with OCD (association p-value < 0.2) were included as covariates. We used the inverse variance method implemented in METAL ( Willer, Li, and Abecasis 2010 ) to meta-analyze summary statistics from each subpopulation and trios for sex-stratified analysis. We performed GWAS and meta-analysis on the combined male/female sample for each subpopulation to ensure that our sex-specific QC yielded results consistent with the recently reported OCD meta-GWAS using the same data ( International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS) 2017 ). The correlation calculated using LD score regression ( B. K. Bulik-Sullivan et al. 2015 ) between our meta-analysis and the previously published meta-analysis ( International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS) 2017 ) was not significantly different from 1 (rg = 1.052, se=0.014). Manhattan and quantile-quantile plots ( Khramtsova and Stranger 2016 ) were used to visualize results. We used PLINK to perform a genotype-sex (GxS) interaction analysis, with principal components as covariates, in each of the individual datasets. We then used METAL to meta-analyze the interaction results. A sex-stratified analysis followed by difference test (Z-score, see below) is equivalent to a formal test for genotype-sex interaction when there is no interaction between covariates and the strata, and the trait variance are equivalent in the two strata. However, different information can be gained from both types of analyses. An interaction test on a combined sample is powered to detect a difference between the sexes in genetic risk and needed to determine whether differences in effect sizes are statistically different between the sexes. On the other hand, a stratified analysis is required in order to characterize the effect size itself, and the direction of effect within each sex. We used Z-scores (correlated with Cochran’s Q statistic but provides directionality of the effect, Supplementary Methods ) to assess heterogeneity between males and females. To obtain a Z-score, and corresponding p-values, for each tested variant, we calculated the differences in effect sizes (beta) between the sexes weighted by the square root of the sum of beta standard errors squared ( equation 1 ). Equation 1 Z − s c o r e = B e t a f e m a l e − B e t a m a l e S E f e m a l e 2 + S E m a l e 2 We define SNPs with Sexually Differentiated Effect (SDEs) as those variants at the extreme ends of the distribution with an absolute value of the Z-score greater than 3 (∣Z-score∣>3), which is roughly equivalent to p<10 −3 , and represents 0.3% of all tested SNPs. We used Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA) SNP2GENE web tool ( Watanabe et al. 2017 ), to perform annotation of the male- and female-specific genome-wide associations. We used the default settings with minor modifications: the minimum p-value of the lead SNP was set at 1.0E-5, the r 2 threshold to define the LD structure of the lead SNPs greater or equal to 0.6, the maximum p-value cutoff at 0.5, MAF>0.01, 250 kb as the maximum distance between LD blocks to merge into a locus, 1000 Genomes Project Phase 3 European population as the reference panel, and variants from the reference panel were included for identification of functional variants in LD with the lead SNP. For the gene-based analysis implemented via generalized-gene set analysis of GWAS data (MAGMA) ( de Leeuw et al. 2015 ) in FUMA, SNPs are mapped to genes if they fall within the gene start and end sites. In this analysis the mean of the chi2-statistic for the SNPs in a gene is calculated and the p-value is obtained from a known approximation of the sampling distribution ( Brown 1975 ; Hou 2005 ). The genome-wide significance threshold is defined as 0.05 / number of genes to which the SNPs are mapped. To calculate the sex-specific narrow-sense SNP-based heritability (h 2 ), (i.e., the proportion of phenotypic variation attributable to the additive effect of all SNP variants in each sex), we used two methods: 1) LD score regression (LDSC) as implemented in LDSC v1.0.0 ( B. K. Bulik-Sullivan et al. 2015 ) and 2) restricted maximum likelihood analysis (REML) implemented in GCTA v1.24.4 ( Yang et al. 2011 ). LDSC analysis was performed on the sex-stratified meta-analysis summary statistics from all study datasets. Meta-analyzed imputed SNPs which overlapped with a panel of high confidence HapMap SNPs were used for the LD score regression. Because our dataset is composed of European individuals, we downloaded precomputed LD scores ( B. K. Bulik-Sullivan et al. 2015 ; B. Bulik-Sullivan et al. 2015 ). Using all individuals, we calculated the total and sex-stratified heritability, checked for residual population stratification (based on the LDSC intercept ( B. K. Bulik-Sullivan et al. 2015 )), and calculated the genetic correlation between males and females. A range of 1–3% OCD population prevalence was used to transform from the observed heritability scale to the liability scale. For REML analysis, we used a combination of the IOCDF-GC and OCGAS European datasets plus the cases from the IOCDF-GC and OCGAS trio dataset and performed an additional PCA analysis on this combined sample to remove any outliers. Genetic relationship matrices (GRM) for autosomes and chromosome X were generated for combined and sex-stratified datasets, removing any individuals who are closely related (IBD>0.05). All pruned imputed SNPs were used to determine the top 20 principal components using smartpca in EIGENSOFT ( Price et al. 2006 ). Genomic-relatedness-based restricted maximum-likelihood (GREML) analysis was performed on the autosomes and the X chromosome (taking into account dosage compensation, Supplementary Methods ) using GRMs and the top 20 ancestry covariates. The same range of prevalence estimates, as in the LDSC analysis, were used to transform heritability to the liability scale. Bivariate GREML analysis was performed to assess the genetic correlation between the sexes. To determine the proportion of the total heritability contributed by each chromosome (including the X chromosome), a separate GRM was generated for each of the 23 chromosomes. Then, all chromosomes were analyzed jointly in a single GREML analysis with 20 PCs to account for population substructure. To assess eQTL enrichment, specifically to test for an enrichment for a gene regulatory role among top GWAS associations and SDEs, we quantified the enrichment of the number of eQTL target genes (eGenes) associated with OCD-associated SNPs. Expression quantitative trait loci (eQTL) enrichment analysis was performed on (a) SNPs nominally associated with OCD (p<10–3) in the combined and sex-stratified GWAS analysis, and (b) SDEs. This analysis is specifically testing the hypothesis that SNPs modestly associated with OCD (within each sex, and within the combined sample) are enriched for immune or brain eQTLs in comparison to null sets of SNPs that are not associated with OCD. Thus, the comparison is only between SNPs associated with OCD and those not associated with OCD (but matched on genomic features). Prior to clumping (r 2 =0.2, 500kb window), each set of SNPs was filtered for variants with fewer than five hundred individuals present in the meta-analysis. We also report results of analyses of unfiltered SNPs ( Supplementary Figure 8 ). eQTL annotation was performed using previously published eQTL results ( Supplementary Table 2 ), including eQTLs derived from 10 regions of the brain and whole blood from GTEx v7 ( GTEx Consortium et al. 2017 ), a meta-eQTL analysis of brain cortex tissue ( Kim et al. 2014 ), as well as CD4+ T cells and CD14+ monocytes ( Raj et al. 2014 ). To assess eQTL enrichment, 1000 randomly ascertained SNPs sets were generated using SNPsnap ( Pers, Timshel, and Hirschhorn 2015 ), sampled without replacement (replacement is allowed only when not enough matched SNPs are available) from the European catalogue of 1000 Genomes SNPs, and matched for minor allele frequency (± 5%), gene density (± 50%), distance to nearest gene (within a 1000kb window), and LD buddies (± 50%) at r 2 =0.8. SNPs in the OCD-associated set and the null matched SNP sets were annotated both with cis -eQTL status and with the genes they regulate (i.e., eGenes) in brain and immune tissues. The enrichment p-value was calculated as the proportion of randomized sets in which the number of eGenes matched or exceeded the observed count among trait-associated SNPs. If multiple variants implicated the same eGene in a tissue or cell type, the eGene was counted only once. This strategy is different from counting individual eQTLs variants, as was done for the previous OCD GWAS ( Stewart et al. 2013 ), where multiple SNPs may be regulating the same gene, while here all eQTLs targeting the same gene are counted only once. We also performed “pan-tissue” eQTL eGene analysis by combining the eQTL results from all the brain tissue subtypes and all the immune tissue and cell subtypes. If an eGene was present in more than one tissue, it was counted only once. To exclude the possibility of eQTL enrichment overestimation due to the gene-rich MHC region, we performed eQTL enrichment analysis both including and excluding SNPs in the HLA region. The enrichment was considered significant if the empirical p-value exceeded Bonferroni multiple testing correction threshold p<0.0036 (i.e. 0.05/14 tissues). We tested for enrichment of anthropometric trait SDEs (ASDEs) among SNPs nominally associated with OCD (p<10 −3 ) in (a) the combined male/female analysis, (b) the sex-stratified analyses, and (c) the OCD SDEs. ASDEs were defined using the approach described in ( Mitra et al. 2016 ) (Z-score p<=10 −3 ) for several anthropomorphic traits from the GIANT consortium ( Randall et al. 2013 ): weight, height, body mass index (BMI), hip circumference (HIP), HIP adjusted for BMI (HIPadjBMI), waist circumference (WC), WC adjusted for BMI (WCadjBMI), waist-to-hip ratio (WHR), and WHR adjusted for BMI (WHRadjBMI) resulting in a total of 12,006 unique ASDEs identified across GIANT phenotypes. We determined the overlap of ASDEs with each OCD subset ( Supplementary Figure 9 ), as well as with 1000 matching SNP sets for each of the OCD subsets. An empirical enrichment p-value was calculated as the proportion of null randomized sets in which the overlap matched or exceeded the observed overlap using the OCD associated SNPs. Genetic correlation analysis of OCD with thirty-one phenotypes of interest was performed for the combined OCD sample and sex-stratified OCD samples using LD score regression ( B. K. Bulik-Sullivan et al. 2015 ). Sex-stratified summary statistics for the following eight phenotypes were obtained ( Supplementary Table 6 ): attention-deficit hyperactivity disorder ( Martin, Walters, Demontis, Mattheisen, Lee, et al. 2017 ), post-traumatic stress disorder ( Duncan et al. 2017 ), reproductive behavior ( Barban et al. 2016 ), insomnia ( Hammerschlag et al. 2017 ), educational attainment ( Okbay, Beauchamp, et al. 2016 ), and alcohol consumption ( Clarke et al. 2017 ; Schumann et al. 2016 ). In the absence of available sex-stratified summary statistics, sex-combined results were obtained for thirty-one phenotypes ( Supplementary Table 7 ): Tourette Syndrome ( Scharf et al. 2013 ; Yu et al. 2015 ), obsessive-compulsive symptoms ( den Braber et al. 2016 ), post-traumatic stress disorder ( Duncan et al. 2017 ), attention deficit hyperactivity disorder ( Neale et al. 2010 ), autism (unpublished, available via Psychiatric Genomics Consortium), bipolar disorder ( Psychiatric GWAS Consortium Bipolar Disorder Working Group 2011 ), major depressive disorder ( Major Depressive Disorder Working Group of the Psychiatric GWAS Consortium et al. 2013 ), schizophrenia ( Schizophrenia Working Group of the Psychiatric Genomics Consortium 2014 ), anxiety disorders ( Otowa et al. 2016 ), depressive symptoms, neuroticism, subjective well-being ( Okbay, Baselmans, et al. 2016 ), anorexia ( Boraska et al. 2014 ), body mass index ( Locke et al. 2015 ), tobacco usage ( Tobacco and Genetics Consortium 2010 ), reproductive behavior ( Barban et al. 2016 ) and structural brain measures (accumbens, amygdala, pallidum, caudate, thalamus, putamen volumes) ( Hibar et al. 2015 ), hippocampal volume ( Hibar et al. 2017 ), intracranial volume ( Adams et al. 2016 ), insomnia ( Hammerschlag et al. 2017 ), educational attainment ( Okbay, Beauchamp, et al. 2016 ), and alcohol consumption ( Clarke et al. 2017 ; Schumann et al. 2016 ). We identified high confidence HapMap SNPs (for which the LD scores have been precomputed) present in the OCD summary statistics and each of the other summary statistics. For continuous traits (e.g. cognitive performance, brain structure volumes) no sample or population prevalence was specified. For binary traits the sample prevalence was calculated based on the reported number of cases in the sample, while the population prevalence was obtained from the literature ( Supplementary Table 7 ).

Genomic control lambda (λ GC ) revealed no significant evidence of population stratification in the male-specific (λ GC =1.019), the female-specific (λ GC =1.026), or the combined (λ GC =1.051) meta-analyses. The intercepts estimated by LD score regression of 1.002, 0.991, 1.005 for sex-combined, female-only, and male-only, respectively, suggested that the mild inflation observed on the quantile-quantile plots (λ GC ) was not due to population stratification but rather to polygenic effects. The Manhattan and quantile-quantile plots ( Figure 1 A-B ) demonstrated no genome-wide significant associations in either males or females. There was little overlap in the top signals across sexes as illustrated by the lack of points along a diagonal line from the left bottom corner to the right upper corner of Figure 1C (top ten associations, Table 1 ). This visualization indicates that there are no variants which strongly associate (-log10(p-value)>3) with OCD in both sexes. In fact, the strongest associations in one sex, have very low -log10(p-value) in the opposite sex. Gene-based tests computed by MAGMA revealed two genome-wide significant genes in the female analysis: GRID2 (pFEMALE = 1.07E-07, pMALE = 7.23E-01) and GPR135 (pFEMALE = 1.55E-06, pMALE = 7.04E-01) which were not significant in males ( Supplementary Figure 4 ). All other tests implemented in FUMA, including MAGMA gene-set analysis, and tissue expression analysis did not result in any significant findings for either sex. The QQ plots ( Supplementary Figure 5 A-B ) of Z-score p-values indicated no significant SDEs (top ten SDEs, Table 2 ), and that the difference in effect size for SDEs was not driven by minor allele frequency (MAF) differences between sexes ( Supplementary Figure 6 ). The MAF distributions for SDEs and all tested SNPs were identical, and sexually-differentiated loci were distributed across the genome proportional to chromosome length ( Supplementary Figure 6D ). P-values from a genotype-sex interaction test ( Supplementary Figure 5 C-D ) were highly correlated with Z-score p-values from the sex-stratified analysis (autosomal SNPs Pearson’s r=0.65, p<2.2e-16, X chromosome SNPs Pearson’s r=0.71, p<2.2e-16). Furthermore, GWAS results in the combined sample with or without sex as a covariate were highly correlated (LDSC r g =0.999, se=0.001). For highly polygenic traits, individual genetic variants, including the most significantly associated variants, typically explain only a small fraction of a trait’s phenotypic variance. To characterize the sex-specific genetic architecture of OCD, we explored sub-threshold associations and their contribution to OCD heritability (h 2 ). The difference in heritability estimates ( Table 3 ) between males (h 2 M =0.131, SE = 0.097) and females (h 2 F =0.296, SE = 0.079), as determined by LDSC regression, was not statistically significant, and the genetic correlation between the sexes was substantial (r g = 1.043, SE = 0.509, p=0.041). The restricted maximum likelihood analysis (REML) estimates of heritability were almost identical between males (h 2 M =0.232, se=0.072, p=0.001) and females (h 2 F =0.240, se=0.057, p=1.07e-05), and to the combined estimate (h 2 =0.238, se=0.033, p=8.621e-14). The REML genetic correlation between males and females was 1.00 (se=0.27). The observed patterns were also robust across population prevalence rates ( Supplementary Table 3 ). One of the mechanisms by which sex differences in OCD could arise is through genetic risk deriving from the sex chromosomes. We observed no significant associations on the X chromosome in either the combined or sex-stratified analyses. A QQ-plot indicated that there was no significant SDEs on the X chromosome ( Supplementary Figure 5 ). Using REML, we estimated the X chromosome (1.6% of total SNPs) heritability as h 2 X =0.010 (se=0.005, p=0.006), which comprised 3.8% of total OCD heritability, and was consistent with expectation ( Supplementary Figure 7 ). When analyzed in each sex separately, X chromosome heritability was not statistically different between females (h 2 FX =0.014, se=0.008, p=0.027) and males (h 2 MX =0.028, se=0.013, p=0.010) at 2.5% OCD prevalence. Results were again robust to estimates derived using a range of OCD prevalence ( Supplementary Table 3 ). To investigate the functional effects of top associations (p<10 −3 ) from the sex-stratified GWAS analysis and SDEs, we annotated each SNP as to whether it was an expression quantitative trait locus (eQTL) in brain or immune tissues. We tested for enrichment of eQTLs derived from brain tissues because brain is the primary tissue of interest, but also eQTLs derived from immune cells because the immune system has been previously implicated in several neuropsychiatric and neurodegenerative traits ( Schizophrenia Working Group of the Psychiatric Genomics Consortium 2014 ; Marsh et al. 2016 ; Heneka, Golenbock, and Latz 2015 ; Furtado and Katzman 2015a ), including OCD ( Furtado and Katzman 2015b ; Murphy, Sajid, and Goodman 2006 ) SDEs showed a significant enrichment for eQTLs from CD4+ T cells (p<0.001), the combination of immune tissues (p<0.002), and combined brain tissues excluding the functionally distinct cerebellum (p=0.001) ( Figure 2 , Supplementary Table 4 ). 188 eGenes were implicated by brain eQTLs and 198 by immune eQTLs, with 48 eGenes deriving from both tissues ( Supplementary Table 5 ). Including HLA SNPs did not significantly affect the results. Consistent with the FUMA results presented above, we did not detect significant enrichments in the combined or sex-stratified analysis. Previous work has revealed enrichment of anthropometric traits SDEs (ASDEs) among top autism (ASD), bipolar disorder (BIP) ( Mitra et al. 2016 ), and endometriosis ( Rahmioglu et al. 2015 ) associated genetic variants, suggesting that the same mechanisms acting on secondary sex characteristic differences later in life may also contribute to sex differences in other complex traits, including neuropsychiatric phenotypes, via pleiotropic effects. There was little overlap and no significant enrichment (p=0.14) for ASDEs among the clumped top combined, female-specific, male-specific GWAS associations, or OCD SDEs ( Supplementary Figure 9 ). As the lower bounds on the genetic correlation estimate of OCD between sexes ranged from 0.49–0.73, we explored whether males and females demonstrate differential genetic correlations between OCD and 30 traits ( Supplementary Table 7 ) which may play a role in OCD development. Our analysis was limited by availability of combined male and female summary statistics for the majority of the traits we tested in the correlation. When available, we have used sex-stratified summary statistics to perform genetic correlation analysis between traits within each sex and assessed differences between sexes. The traits chosen for analysis included (1) neuropsychiatric phenotypes and behavioral traits (many of which exhibit sexually-differentiated characteristics), (2) traits which overlap with known sexually-differentiated clinical symptoms in OCD (e.g. smoking, eating disorders-anorexia, and body mass index), (3) brain structure volumes, and (4) reproductive behavior (age at first birth and number of children ever born). Using sex-stratified summary statistics for OCD and traits listed in Supplementary Table 6 , we performed cross-trait genetic correlations within each sex (i.e. Rg between male OCD with male ADHD and female OCD with female ADHD). There were no significant cross-trait genetic correlations in either sex after multiple testing correction nor did the cross-trait correlations differ significantly between males and females ( Supplementary Table 6 ), suggesting that sex stratified analyses are likely still underpowered. To increase the power of this analysis, next, we used sex-combined summary statistics for several complex traits ( Supplementary Table 7 ) and performed cross-trait correlation analysis with sex-stratified OCD summary statistics. Several traits (bipolar disorder, schizophrenia, and neuroticism) exhibited a significant genetic correlation with female OCD, but not male OCD, again, possibly influenced by sample size ( Supplementary Table 7 ). Again, the genetic correlations for OCD with other traits did not differ significantly between males and females. A genetic correlation analysis using sex-combined summary statistics for OCD and other traits, revealed several significant cross-trait correlations, indicating that a larger sample size results in more precise estimates of the genetic correlation. We observed novel significant genetic correlations between the sex-combined OCD sample and the sex-combined summary statistics from age at first birth (r g =0.37, se=0.07, 4.83e-07), number of children ever born (r g =−0.35, se=0.09, p=6.66e-05), and replicated previously published results ( Brainstorm Consortium et al. 2018 ; Davis et al. 2013 ; Yu et al. 2015 ).

Obsessive-compulsive disorder is one of many neuropsychiatric traits exhibiting sex differences in both age of onset and presentation of symptoms. Gene-based analysis identified two genes (GRID2 and GPR135) with female-specific associations that were not present in males, however, at the level of individual loci, no genome-wide significant associations were detected in either the sex-stratified GWAS or the genotype-sex interaction analysis. The genome-wide genetic correlation for OCD between males and females was not significantly different from 1 and OCD heritability estimates were not significantly different between the sexes. Additionally, we observed no significant differences in the cross-trait genetic correlations between males and females which is currently best explained by the absence of ubiquitous genetic architecture differences between male and female OCD, as well as small sample sizes which negatively impact on the ability to detect smaller differences between the sexes. Finally, partitioned heritability analysis indicated that the X chromosome contributed to the polygenic liability of OCD, underscoring the importance of including the X chromosome in GWAS of OCD. The GRID2 gene is part of the glutamatergic signaling system ( Pittenger, Bloch, and Williams 2011 ) which is thought to be important in OCD and is expressed in the brain regions which have been implicated in OCD (cerebellum, caudate, putamen, nucleus accumbens, and the anterior cingulate cortex) ( Graybiel and Rauch 2000 ). Less is known about GPR135, however, results from the GTEx portal ( https://gtexportal.org ) indicate that it is also expressed in brain. Taken together, these results indicate that significant sex-specific effects for OCD likely exist but will be challenging to detect given their modest effect sizes and the sample size required to detect statistically robust genotype-sex interactions. Encouragingly, significant sex-stratified associations have been identified in studies of ASD and ADHD, demonstrating the value of increasing sample size for the study of sexually-differentiated genetic effects ( Mitra et al. 2016 ; Martin, Walters, Demontis, Mattheisen, Lee, et al. 2017 ). Furthermore, we observed that SDEs, SNPs with the greatest heterogeneity in effect size between males and females were enriched for gene regulatory function (eQTLs) in brain and immune tissues, implicating these tissues in sexual-differentiation of OCD. The enrichment of immune eQTLs among SDEs is consistent with both the known role of the immune system in several neuropsychiatric traits ( Schizophrenia Working Group of the Psychiatric Genomics Consortium 2014 ; Furtado and Katzman 2015a , [b] 2015 ; Murphy, Sajid, and Goodman 2006 ), and the observed sex differences in immune function ( Klein and Flanagan 2016 ). Moreover, a recent whole-exome sequencing study found that OCD probands have a higher rate of de novo nonsynonymous single-nucleotide variants in genes enriched for neurodevelopmental and immunological processes ( Cappi et al. 2016 ). Though specific mechanisms remain unknown, these studies provide new evidence for an old hypothesis linking the immune system with compulsive behavior and OCD ( Marazziti et al. 1999 ; Kawikova et al. 2007 ; Murphy, Sajid, and Goodman 2006 ; Slattery et al. 2004 ; Miguel et al. 1995 ; Swedo et al. 1998 ; Murphy et al. 2012 ; Snider and Swedo 2004 ; Swedo et al. 2012 ; Leonard et al. 1992 ; Carapetis and Currie 1999 ). Limitations for this study include sample size, and ascertainment strategies that may bias towards earlier age of onset which could result in uneven representation of disease subclasses among males and females. For example, early-onset OCD is reportedly slightly more heritable than adult-onset ( Davis et al. 2013 ; G. Nestadt et al. 2000 ; van Grootheest et al. 2005 ). Thus, uneven representation of males and females in the early- and adult-onset OCD groups could confound heritability if estimates are influenced by both sex and age-at-onset. Age-of-onset information is incomplete in many of the historical sample collections that have been included in this meta-analysis. The lack of detailed clinical data limits our ability to address many important questions related to symptom type, symptom severity, and age of onset. These limitations underscore the need for larger OCD datasets phenotyped in greater detail to delve deeper into both genetic and clinical sex differences observed in OCD.

Obsessive-compulsive disorder (OCD) displays sex differences in age of onset, progression, and symptomatology, however, the genetic basis of sex differences in OCD has not yet been comprehensively explored ( Flament et al. 1990 ; Swedo et al. 1989 ; Bellodi et al. 1992 ; Boileau 2011 ). Epidemiological studies indicate a worldwide lifetime prevalence of OCD between 1 and 3% ( Kessler et al. 2005 ; Ruscio et al. 2010 ; Torres and Lima 2005 ; Weissman et al. 1994 ) and while boys comprise approximately two thirds of the childhood cases of OCD, typically defined as onset before age 15 ( Flament et al. 1990 ; Swedo et al. 1989 ; Bellodi et al. 1992 ; Boileau 2011 ), females predominate the late-onset cases of OCD. In addition to demonstrating a later age of onset, females with OCD have higher rates of precipitating events which include pregnancy and childbirth ( Lochner et al. 2004/3 ; Gerald Nestadt, Grados, and Samuels 2010 ). Compared to females, males with OCD report more religious, sexual, and symmetry symptoms, more alcohol dependence, and lower rates of marriage and employment. Females with OCD are more likely to be married, report more sexual abuse during childhood, often report exacerbation of symptoms in the premenstrual/menstrual period, during/shortly after pregnancy, with menopause, and tend to have more contamination and cleaning compulsions, as well as eating disorders, reviewed in Mathis et al ( Mathis et al. 2011 ). Although the recently-published genome-wide association studies of OCD ( Mattheisen et al. 2015 ; Stewart et al. 2013 ; International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS) 2017 ) found no genome-wide significant associations, these studies demonstrated that common variants account for a significant proportion of OCD heritability (24–32%) ( Davis et al. 2013 ; International Obsessive Compulsive Disorder Foundation Genetics Collaborative (IOCDF-GC) and OCD Collaborative Genetics Association Studies (OCGAS) 2017 ). They also indicate that the strongest associated variants in OCD genome-wide association study (GWAS) are enriched for expression quantitative trait loci (eQTLs) and methylation QTLs derived from frontal lobe, cerebellum, and parietal lobe tissue ( Stewart et al. 2013 ), demonstrating that biologically meaningful associations exist within the top ranked SNPs and that increasing sample sizes will likely identify significant common variant associations for OCD risk. In addition to increasing sample size, another approach to improve power for GWAS is to reduce genetic heterogeneity. For example, if sex significantly modifies the effect of genetic variation on the risk of OCD, then combining males and females with OCD may weaken or obscure sex-specific effects. For example, previous studies have discovered novel loci which were previously undetected due to heterogeneity between sexes ( Mitra et al. 2016 ; Martin, Walters, Demontis, Mattheisen, Hong Lee, et al. 2017 ; Taylor et al. 2013 ; Randall et al. 2013 ; Liu et al. 2012 ; Winkler et al. 2015 ; Hartiala et al. 2016 ; Orozco et al. 2012 ; Zhuang and Morris 2009 ; Singh et al. 2016 ). Motivated by the sex differences in OCD, we tested the hypothesis that the genetic architecture of OCD varies between the sexes. We first performed a sex-stratified genome-wide association meta-analysis and genotype-sex interaction meta-analysis including autosomes and the X chromosome. We then developed an approach to identify SNPs with Sexually Differentiated Effects (SDEs), and assessed whether the SDEs regulate gene expression and are enriched for associations with sexually-differentiated anthropometric traits (i.e. height, weight, body mass index, hip and waist circumference) as observed in autism spectrum disorders ( Mitra et al. 2016 ). Third, we performed SNP-based heritability analysis to (a) assess the proportion of overall OCD heritability explained by the X chromosome, and (b) test for evidence of sexually-dependent liability threshold for OCD between males and females. An important correlate of the sex-dependent liability threshold is that the sex with the lower prevalence/milder presentation requires a higher genetic burden to become affected and therefore is more likely to have affected children, also known as the Carter effect. This has been reported for several complex traits ( Kruse et al. 2012 ; Kantarci et al. 2006 ). Fourth, we performed a sex-stratified genetic correlation analysis with other traits which may play a role in OCD development (e.g. brain volumes), are sexually-differentiated (e.g. autism, Tourette syndrome, attention deficit hyperactivity disorder, etc.), or are known to show differences in comorbidity between males and females with OCD (e.g. smoking, eating disorders, and reproductive behavior). Here, we present the first genome-wide assessment of the sex-specific genetic architecture of OCD utilizing the largest OCD dataset currently available. We also provide best practices for sex-stratified analysis which can be adopted in future studies of OCD and other phenotypes.