{"paper_id":"249251f4-5ac8-4c1a-ae5e-f58ee27fc333","body_text":"Genetic effects on transcriptional regulation underlie the pathogenic mechanisms of many human traits and diseases ( Peters  et al ., 2016 ;  Ongen  et al ., 2017 ;  Gamazon  et al ., 2018 ). Genetic variants that regulate gene expression, termed expression quantitative trait loci (eQTLs) can be shared between tissues or can be tissue specific ( Consortium  et al ., 2017 ). Consequently, to investigate genetic contributions to disease mechanisms, we must investigate the tissue, and ultimately individual cell types, relative to a disease to understand the genetic contribution to pathogenesis.\nThe endometrium, as the inner most layer of the uterus is a vitally important reproductive tissue, integral for fertility and implicated in many reproductive disorders. It undergoes cyclical changes largely driven by hormonal regulation and changes in cellular composition ( Evans  et al ., 2016 ). Our previous studies have shown endometrial gene expression and methylation is influenced by genetic variation ( Powell  et al ., 2016 ;  Fung  et al ., 2017 ;  Fung  et al ., 2018 ;  Mortlock  et al ., 2019 ) and suggest these genetic variants contribute to the susceptibility to reproductive disorders. The specific genes that are influenced by these genetic variants and that contribute to disease pathogenesis however remain to be elucidated.\nEndometriosis, characterised by endometrial-like tissue that form lesions outside the uterus is a common reproductive disorder affecting 6–10% of reproductive aged women and it is believed to stem from endometrial tissue ( Bulletti  et al ., 2010 ;  Giudice, 2010 ). A recent genome-wide association study (GWAS) identified 27 genomic loci associated with endometriosis ( Rahmioglu  et al ., 2018 ) and our previous studies showed that expression of critical endometrial target genes and methylation of CpG sites are altered in genomic regions associated with endometriosis susceptibility ( Powell  et al ., 2016 ;  Fung  et al ., 2017 ;  Fung  et al ., 2018 ;  Mortlock  et al ., 2019 ). This suggests that susceptibility to endometriosis is mediated by changes in endometrial gene expression and methylation under the control of genetic risk factors. The underlying mechanisms increasing endometriosis susceptivity at many of these genetic regions is not yet clear.\nRNA sequencing (RNA-seq) is a powerful gene expression technique that quantitates individual RNA transcripts with a broader dynamic range than microarray technology, which captures only 30% of the data available in RNA-seq ( Mortazavi  et al ., 2008 ;  Wang  et al ., 2009 ;  Zhao  et al ., 2014 ). The aim of this study was therefore to extend our understanding of the genetic regulation of transcription in endometrium by analysing paired-end total RNA sequence data from 206 endometrial samples, characterize the eQTLs in the endometrium and their similarity to other tissues, and determine their association with susceptibility genes for reproductive traits and diseases, such as endometriosis.\n\nTo identify the tissue types that could be associated with endometriosis, we performed cell-type enrichment analysis using the method outlined in  Finucane  et al . (2018) . We used the summary statistics from the endometriosis meta-analysis conducted by  Sapkota  et al . (2017) , which contained 17 045 endometriosis cases and 191 596 controls. Gene expression data were obtained from the GTEx project ( Consortium, 2015 ,  Consortium  et al ., 2017 ) and the Franke Lab ( Pers  et al ., 2015 ). The GTEx dataset contains gene expression data from RNA-Seq analysis of 53 different tissues or cell types from 8550 human samples. The Franke Lab gene expression data comes from microarray analysis of 152 distinct tissues and cells types from 37 427 human samples. The different tissues or cell types were classified into nine groups for the visualisation, including adipose, blood or immune, cardiovascular, central nervous system, digestive, liver, musculoskeletal-connective, pancreas and other.\nNumber of patients recruited from the RWH and Melbourne IVF Clinic (IVF) and associated endometriosis and stage of menstrual cycle information.\nFor each gene in all of the distinct tissue or cell-type gene expression datasets, we calculated the  t -statistic, a measure of tissue specific gene expression. Within each tissue or cell type, the  t -statistics for each gene was ranked and only the top 10% were used in the next step of the analysis. From each gene, we added a 100-kb window on either side of the transcribed region to construct the genome annotation and then performed stratified linkage disequilibrium (LD) score regression on the endometriosis GWAS summary statistics to test whether the disease heritability was enriched in loci containing genes with the highest expression in particular tissues. Tissue or cell types with FDR < 0.05 were classified as significantly enriched.\nWoman of European ancestry and reproductive age were recruited from clinics at the Royal Women’s Hospital (RWH) and Melbourne IVF Clinic (IVF) in Melbourne, Australia. A total of 206 women were included in the study, consisting of 184 RWH gynaecology patients and 22 IVF patients. RWH patients underwent investigative laparoscopic surgery in response to pathological symptoms and/or infertility during which endometrial tissue samples were extracted by curettage. A detailed clinical history and surgical and pathological results were obtained for each participant in RWH sample set.\nThe endometriosis status was recorded following surgical diagnosis at laparoscopy for women from RWH, or self-reported for women from the IVF clinic. A histological assessment of each endometrial biopsy was performed by an experienced pathologist to categorise samples into seven different menstrual cycle stages including menstrual (M), early-proliferative (EP), mid-proliferative (MP), late-proliferative (LP), early-secretory (ES), mid-secretory (MS) and late-secretory (LS). We excluded samples that were from women of non-European ancestry or who underwent hormonal treatment and tissues that showed abnormality on histopathology or with ambiguous diseases status or cycle stage. A summary of samples and associated clinical detail is provided in  Table I . From each individual, we collected whole blood samples prior to surgery and endometrial tissue. The endometrial tissue was stored in RNAlater (Life Technologies, USA) at −80°C for later RNA extraction, while whole blood was prepared for DNA isolation.\nThe study was approved by the Royal Women’s Hospital Human Research Ethics Committee (Projects 11–24 and 16–43), the Melbourne IVF Human Research Ethics Committee (Project 05–11) and the University of Queensland. Informed consent was obtained from all participants.\nDNA samples from each of the 206 individuals were genotyped on HumanCoreExome or Infinium PsychArray chips (Illumina, USA). Quality control (QC) was performed in PLINK according to the protocol outlined in  Fung  et al . (2018) . Following QC, a total of 282 625 SNPs (hg19) were phased using Shapelt V2 and taken forward to imputation using the haplotype reference consortium reference panel (version r1.1 2016) on the Michigan Imputation Server. SNPs with low imputation quality (R 2  < 0.8), missing rate > 5%, minor allele frequency (MAF) < 1 × 10 –4,  and Hardy–Weinberg equilibrium < 1 × 10 −6  after imputation were removed. The remaining SNP positions were lifted-over to the Ensembl genome build 38 (GRCh38) using CrossMap v.0.2.8. SNPs failing to lift-over were assigned to their new GRCh38 position manually based on dbSNP151 GRCh38 patch release 7 (GRCh38.p7), leaving 6 230 993 SNPs for further analysis.\nTotal RNA was isolated from endometrial tissue using the Allprep DNA/RNA Mini Kit (Qiagen, CA) as per the manufacturer’s instructions. RNA quality was checked using the Bioanalyzer 2100 (Agilent Technologies, CA) and RNA concentration was measured using the NanoDropND-6000 (Thermo Fisher Scientific, USA). All samples were high quality with an RNA integrity number greater than 8.\nThe RNA samples were treated with Turbo DNA-free kit (Thermo Fisher Scientific, USA) prior to RNA-seq library generation. Stranded RNA-seq libraries were prepared using the Illumina TruSeq Stranded Total RNA Gold protocol which includes ribosomal depletion (Illumina, USA). Libraries were pooled and sequenced to a mean depth of 37 490 673 reads (178 samples; 75 bp pair-ended reads) on the Illumina HiSeq 4000 and 40 818 062 reads (28 samples; 120 bp pair-ended reads) on the Illumina HiSeq 2000 (Ilumina, USA). Raw sequencing reads were quality checked using FastQC v0.11.7( Andrews, 2010 ) and MultiQC v1.6( Ewels  et al ., 2016 ). Low quality reads and contaminating HiSeq Illumina adapter sequences were trimmed using Trimmomatic v0.36( Bolger  et al ., 2014 ).\nTrimmed reads were aligned against the human reference genome (Ensembl  Homo sapiens  GRCh38 release 84) using HISAT2 v2.0.5. Transcript assembly was performed using StringTie v1.3.1( Pertea  et al ., 2015 ,  2016 ) and the Ensembl  Homo sapiens  GRCh8 release 91 reference assembly. Reads mapping to each known transcript were directly counted in StringTie to generate transcript-, exon- and intron-level expression matrices in ‘fragments per kilobase of transcript per million mapped reads’ units for each individual. Raw gene count matrices were also produced using a Python script provided by StringTie.\nGenes expressed at a low level, i.e. genes with counts per million (CPM) <0.22 (~10 counts) and expressed in <90% of the samples, were removed. Raw gene counts were normalized for composition bias and total raw reads (library size) using the Trimmed Mean of M (TMM) ( Fadista  et al ., 2014 ;  Taneera  et al ., 2015 ;  Seo  et al ., 2016 ) method in the edgeR R package v3.22.3 ( Robinson  et al ., 2010 ). Normalized counts were converted to CPM and log2 transformed (log2-CPM).\nRNA-seq counts and genotype data from the 206 individuals was used to test association between genotype and gene expression. A total of 17 022 genes expressed in >90% samples were included in the analysis. Expression values were TMM-normalized, converted and log2 transformed (log2-CPM) as described above. Individual level genotype data for 6 230 993 SNPs was also included in the analysis. The  cis -eQTL analysis was carried out using a linear regression model in the MatrixeQTL R package v2.2( Shabalin, 2012 ).  cis- eQTLs were defined as SNPs located within ±250 kb from gene start and stop position. Batch effects (lanes within the same flow-cell and between flow-cells), stage of menstrual cycle and endometriosis diseases status where included in the model as covariates. Comparisons with SNPs on a different chromosome to the associated gene, were classified as  trans -eQTLs. The  trans -eQTL analysis was also performed using MatrixeQTL with the same covariates, this time setting a MAF threshold of >0.05 leaving 4 922 014 SNPs to be included in the analysis. To identify independent  cis -eQTL signals, we performed a conditional analysis using the COJO method in GCTA ( Yang  et al ., 2011 ;  Yang  et al ., 2012 ), including the effect of the eSNP (SNP associated with gene expression) with the smallest  P -value as a covariate in the model for each of the 444 genes with Bonferroni significant eQTLs.\nTo evaluate whether the genetic effects on gene expression in endometrial tissue also occurred in other tissues, we correlated eQTL effects with 48 tissues from the GTEx v7 project ( Supplementary Table SI ) ( Consortium, 2015 ;  Consortium  et al ., 2017 ) and the blood eQTL dataset from eQTLGen Consortium (eQTLGen) consisting of 31 684 individuals ( Võsa  et al ., 2018 ). We used the  r b  method developed by  Qi  et al . (2018)  to estimate the correlation between genetic effects at top  cis -eQTLs whilst accounting for eQTL effect estimation errors ( Qi  et al ., 2018 ). Briefly, we used the top significant brain eQTLs ( P eQTL  < 5 × 10 −8 ) from the religious orders study and memory and ageing project (ROSMAP) ( Ng  et al ., 2017 ) as a reference to avoid ascertainment bias. Subsequently top ROSMAP  cis -eQTLs present in both the endometrium eQTL dataset and the GTEx tissue eQTL dataset being tested were used in the effect size (ES) correlation analysis. Genome positions of endometrial eQTLs from this study were converted from the GRCh38 assembly back to the GRCh37 assembly for direct comparison with ROSMAP, GTEx and eQTLGen datasets.\nTo identify potential endometrial specific  cis -eQTLs, we examined the overlap of Bonferroni significant, independent endometrial  cis -eQTLs with  cis -eQTLs in the 48 different tissue types in the GTEx ( Supplementary Table SI ) and eQTLGen datasets ( Võsa  et al ., 2018 ).  cis -eQTLs were matched between tissues based on the same eSNP and gene associations including eSNPs in linkage with the lead eSNP (r 2  > 0.8). The direction of effect was also considered for those with the same eSNP.\nThe limma R package was used to analyse genome-wide differential gene expression following the removal of genes expressed at a low level. Only genes with a minimum of 10 counts and expressed in at least 90% of samples were analysed for differential expression. We performed three distinct differential gene expressions comparisons: (i) between endometriosis cases and controls, (ii) between cases and controls in the mid-proliferative stage and (iii) between cases and controls in the mid-secretory stage. We fitted batch effect (flow-cell and lane) as covariates in all three models and corrected for stage of menstrual cycle in model 1. The normalized counts were transformed using the voom ( Law  et al ., 2014 ) function in limma before being fitted to the linear model. The eBayes method was used to contrast between groups. Resulting  P -values were adjusted using Benjamini–Hochberg with a significance threshold of 0.05. Only genes passing FDR < 0.05 were categorized as being significantly differentially expressed genes.\nWe assessed the overlap of endometrial eQTLs and endometriosis associated SNPs. To date, 27 loci have been associated with endometriosis based on a GWAS meta-analysis across different populations ( Sapkota  et al ., 2017 ;  Rahmioglu  et al ., 2018 ). To identify overlaps between endometrial  cis- eQTLs and endometriosis risk loci, we analysed LD between significant  cis- eQTL (FDR < 0.05) eSNPs and the lead risk SNPs within the 27 risk loci. eQTLs with r 2  > 0.8, as calculated using LDlink ( Machiela and Chanock, 2015 ), were defined as being in LD with the related risk SNP. We also assessed the overlap with SNPs associated with various other traits in the GWAS catalogue.\nIn the absence of gene expression data from large cohorts, new powerful statistical approaches allow us to impute gene expression and conduct a transcriptome-wide association study (TWAS). The resulting association statistics reflect underlying relationships between gene expression and disease risk contributing to the identification of target genes underlying complex traits. Transcriptome-integrated genetic association resource (TIGAR) is a software tool used to impute transcriptomic data and perform TWAS using summary-level GWAS data. We used the Train Dirichlet process regression imputation model in TIGAR to estimate  cis -eQTL ES in endometrium using the 206 RNA-seq samples with matched imputed genotype data. Using estimates from the training model as weights, we conducted a TWAS using endometriosis GWAS summary-statistics from the  Sapkota  et al.  (2017)  meta-analysis. Details of the training and association models can be found in  Nagpal  et al . (2019) . To investigate multiple genes in significantly associated loci, we tested for a correlation between the predicted gene expression in each loci. Gene expression was imputed using estimated effects sizes from the training model and individual genotype data from 5186 individuals belonging to the QIMRHCS cohort ( Nyholt  et al ., 2012 ;  Sapkota  et al ., 2017 ). Gene expression for each individual was predicted by matrix multiplication between estimated ES and additive genotype (X g ) (Predicted_GE = ES * X g) . The Pearson correlation was calculated between the predicted expression of each gene pair in a locus using the cor function in the R environment.\nTo assess the association between genetic variant, gene expression and trait, we performed summary data-based Mendelian randomization (SMR) ( Zhu  et al ., 2016 ). The SMR was conducted by integrating the summary eQTL data and the GWAS meta-analysis from  Sapkota  et al . (2017) . Associations passing the SMR test and with a  P HEIDI  of >0.05/(number of genes passing the SMR test) were considered significant. Due to the limited power of our eQTL analysis, the presence of both secondary eQTL signals and secondary signals in multiple endometriosis risk loci, we also conducted a multi-SNP-based SMR to use information from all significant SNPs in each region ( Wu  et al ., 2018 ). To further avoid confounding effects of multiple signals, reduce multiple testing burden and leverage the power of the large GWAS, we conducted a modified SMR analysis in which the trait (endometriosis) was treated as the exposure and the gene expression within each GWAS locus (2 MB either side of top SNP) was treated as the outcome. Using this modified SMR approach the top GWAS SNPs are selected for testing based on their association with endometriosis. To increase power, we repeated the SMR analyses using the large eQTLGen blood eQTL dataset as a proxy ( Võsa  et al ., 2018 ). Using available GWAS summary statistics for age at menopause ( Day  et al ., 2015 ) and epithelial ovarian cancer ( Phelan  et al ., 2017 ), we conducted an SMR analysis using endometrial eQTLs against these traits.\n\nUnderstanding the role of genetic regulation of gene expression in complex traits and diseases is more powerful when measuring genetic effects in disease relevant tissues. Using summary statistics from  Sapkota  et al . (2017)  and gene expression data from a range of tissues and cell types, we show a significant enrichment for genes annotated to endometriosis-associated loci in female reproductive tissues, myometrium, arteries, serum and digestive tissue ( Table II ,  Fig. 1 ). This suggests effects of regulation of gene expression on endometriosis risk loci are more likely to be detected in reproductive tissues including the endometrium.\nNumbers preceding tissues correspond to the National Institute of Health (NIH) medical subject heading (MeSH) tree structure numbers used to label and distinguish tissues in the Franke lab dataset.\nFollowing assembly and quantification of RNA-seq reads in endometrial samples from 206 European women, we identified 29 791 genes expressed (CPM > 0.22; ~10 counts) in at least two samples. On average samples expressed 25 980 genes. When restricting the gene set to only genes expressed in >90% of samples, we retained 17 022 Ensembl genes. The 12 769 genes only expressed in 1–90% of samples were expressed in varying proportions of samples reflecting the complex structure of gene expression in endometrium previously reported ( Fung  et al ., 2018 ).\nIntegrating data for 6 230 993 (4 922 014 for  trans ) genotyped and imputed SNPs with the RNA-seq data from the 206 European women (17 022 genes expressed in >90% samples), and following Bonferroni correction for multiple testing, we detected genetic effects on gene expression in endometrium for 444 sentinel  cis -eQTLs ( P  < 2.57 × 10 −9 ) and significant  trans -eQTLs for 30 genes ( P  < 5.97 × 10 −13 ) ( Table III ) ( Supplementary Fig. S1a and b ). An additional 22 independent secondary  cis -eQTL signals were detected following conditional analysis whereby the association test was rerun for each gene conditioning on the effect of the primary sentinel SNP. We identified novel  cis -eQTLs, following Bonferroni correction, for 327 genes not previously reported as significant in endometrium. The eQTL dataset generated in this study is available at  http://reproductivegenomics.com.au/shiny/endo_eqtl_rna/ .\nWhen comparing these results to our previous endometrium microarray study ( Fung  et al ., 2018 ), we find 75% of Bonferroni significant eQTLs identified by RNA-Seq were nominally significant ( P  < 0.05) in the microarray data with ES highly correlated (R = 0.75). This reflects the high correlation between the average expression of each gene in endometrial tissue measured by RNA-Seq versus microarray (R = 0.902 ± 0.01) ( Supplementary Fig. S2a–e ) and a high correlation of gene expression within samples (R = 0.789 ± 0.101) ( Supplementary Fig. S3 ). There were 28 FDR significant  trans -eGenes (genes with a trans-eQTL) replicated between the two endometrial datasets.\nOverall, we observed good global correlation in genetic effects on gene expression between tissues, suggesting that many eQTLs exhibit more general shared effects. The correlation in eQTL effects ( r b ) between endometrium and the 48 tested tissues ranged from  r b  = 0.54 with brain cerebellum to  r b  = 0.72 with minor salivary gland ( Fig. 2a ,  Supplementary Table SII ), where  r b  represents the correlation of eQTL effects between tissues accounting for estimation errors in the eQTL effects ( Qi  et al ., 2018 ). Female reproductive tissues including vagina, uterus, breast mammary tissue and ovary all had correlations above 0.68 ( Fig. 2a and b ). We observed large correlations between genetic effects in endometrium and digestive tissues such as the salivary gland ( r b  = 0.72), stomach ( r b  = 0.72) and colon ( r b  = 0.71) and individual cell types such as fibroblasts ( r b  = 0.72) ( Fig. 2a and b ). A lower correlation was observed with central nervous system tissues and hematopoietic tissues ( Fig. 2a and b ). The genetic effects on gene expression in endometrium had the highest correlation with tissue derived from the endoderm (average  r b  = 0.69) ( Supplementary Fig. S4 ), compared to tissue derived from the other germ cell layers.\n* The number of  cis -eQTLs and  trans -eQTLs that pass the less stringent Benjamini-Hochberg threshold of FDR <0.05 is included for comparison (P cis  < 5.32 × 10 −4 ; P trans  < 4.73 × 10 −8 )\nA total of 305 Bonferroni significant endometrial eQTLs were also reported in at least one other GTEx tissue, with an average of 18 reported in any single tissue. The number of endometrial eQTLs reported in each GTEx tissue was highly correlated with the sample size (correlation = 0.88). When comparing endometrial eQTLs with the 48 GTEx tissues, we found a large proportion showed high tissue specificity, observed in only three tissues or fewer (174 sentinel eQTLs). Alternatively, the eQTLs found in more than three tissues were commonly observed across most tissues (185 sentinel eQTLs in >24 tissues) ( Supplementary Fig. S5 ). However, when checking for overlap with the much larger eQTLGen blood dataset ( n  = 31 684) ( Võsa  et al ., 2018 ), 71.6% of endometrial eQTLs were also reported in blood. The large overlap between eQTLs identified in endometrium and blood is likely due to the large sample size in the eQTLGen dataset which has sufficient power to detect eQTLs with much smaller effects. We identified 68 sentinel  cis -eQTLs that may be specific to endometrium ( Supplementary Table SIII ), although a large proportion (85%) of the  cis -eQTLs in endometrium were also reported in other tissues.\nDifferential gene expression between endometriosis cases ( n  = 143) and controls ( n  = 63) was tested using limma in R. Following correction for stage of cycle, no genes were significantly differentially expressed. Differences between cases and controls within proliferative and secretory stages with the largest sample sizes was also tested. No significant differences were observed between cases ( n  = 56) and controls ( n  = 16) in the mid-proliferative stage. Nominal differences in expression between cases ( n  = 27) and controls ( n  = 14) were observed for 43 genes in the mid-secretory stage ( Table IV ). Expression differences for these 43 genes however, were no longer significant following correction for the multiple testing of the three different comparisons, suggesting that these may represent false positives or that there is insufficient power to detect subtle genome-wide significant effects in our current dataset.\nGenes differentially expressed between women with and without endometriosis in endometrium in the mid-secretory stage of the menstrual cycle. Genes with a positive fold change are upregulated in cases and those with a negative fold change are downregulated in cases.\nWe performed a transcriptome-wide association analysis to identify gene expression associated with endometriosis risk. Gene expression and genotype data from the 206 samples was used to estimate the weighted effect of each SNP on each  cis -gene and combined with summary-level endometriosis GWAS data ( Sapkota  et al ., 2017 ) to impute gene expression and perform a TWAS. Using a transcriptome-wide significance threshold of 3.28 × 10 −6 , we identified 252 genes associated with endometriosis located at 39 independent loci ( Fig. 3 ,  Table V ,  Supplementary Table SIV ). Five of these loci harboured genome-wide significant SNPs associated with endometriosis including; rs1903068 on chromosome 4 near kinase insert domain receptor ( KDR) , rs12700667 on chromosome 7 near  LOC100506236 , rs10090060 on chromosome 8 near ganglioside induced differentiation associated protein 1 ( GDAP1 ), rs1802669 on chromosome 10 near MLLT10 Histone Lysine Methyltransferase DOT1L Cofactor ( MLLT10 ) and rs4762326 on chromosome 12 near vezatin ( VEZT ) ( Table V ). The remaining loci all included nominally significant SNPs from the GWAS ranging from  P  = 2.64 × 10 −3 –4.40 × 10 −7  ( Table V ). Many implicated regions contain multiple significant genes whose predicted gene expression was, in most cases, correlated ( Supplementary Fig. S6 ). The most significant gene at each locus is presented in  Table V . This analysis highlights genes potentially involved in endometriosis pathogenesis previously associated with endometriosis risk and provides stronger support for novel regions showing only nominal significance in GWAS results.\nGWAS data is from  Sapkota  et al . (2017) . All positions are based on the hg19 genome version.\nLoci associated with endometriosis and the most significant GWAS SNPs within these regions.\nFunctional annotation of genetic variants associated with complex traits and diseases is another valuable method to identify target genes and prioritise them for further study. Previous studies have reported eQTLs associated with endometriosis on chromosome 1 ( Powell  et al ., 2016 ) and 12 ( Holdsworth-Carson  et al ., 2016 ). In this study, we identified eQTLs for two genes,  VEZT  and FYVE, RhoGEF and PH domain containing 6 ( FGD6 ) overlapping the GWAS signal for endometriosis risk locus on chromosome 12 ( Table VI ). SMR is a method for testing whether the eQTLs and GWAS signals overlap by chance or have some causal association. Implementation of a standard SMR analysis ( Zhu  et al ., 2016 ) found no significant associations between genetic variants, endometrial gene expression and endometriosis risk. The standard implementation for SMR analysis tests only the most significant eQTL at each locus and assumes only one causal variant is associated with both gene expression and the trait. To assess the possibility that additional independent eQTLs and GWAS signals in each region contribute to heterogeneity and dilute the signal, we conducted a multi-SNP-based SMR (SMR-multi) ( Wu  et al ., 2018 ). This increases the power to detect a pleiotropic signal by including multiple SNPs within each individual  cis -eQTL locus. There were no signals that passed the SMR-multi analysis.\nA modified SMR selecting the most significant endometriosis-associated GWAS SNPs as the instrument, thereby treating endometriosis as the exposure, reduced the multiple testing burden and detected three signals that pass both the SMR test and the HEIDI (HEterogeneity In Depedent Instruments) test;  FGD6, VEZT  ( Fig. 4 ;  Table VII ) and  AL022068.1  ( Fig. 5 ;  Table VII ). The HEIDI test is used to distinguish independent overlapping signals and the same SNP influencing expression and disease risk. SNPs with low HEIDI  P -values have a higher probably of being independent overlapping signals and are subsequently disregarded. The expression of both  FGD6  and  VEZT  was also significantly associated with endometriosis in the TWAS analysis.\nP-value significance thresholds:\nP SMR  < 9.62E-04\nP HEIDI  > 0.017\nResults of modified SMR analysis using endometrial eQTLs and endometriosis summary statistics.\nTo increase power, we performed the same SMR analysis using the larger eQTLGen blood eQTL dataset as a proxy. Whilst this increases the power to detect genetic variants that effect gene expression consistently in blood and endometrium, we lose power to detect endometrium specific signals. Using the standard SMR approach, we detected one association between rs2473290 and the cell division cycle 42 ( CDC42 ) gene that passed the SMR test ( P SMR  = 5.77 × 10 −10 ), but did not pass the HEIDI test ( P HEIDI  = 3.74 × 10 −5 ). We repeated the SMR analysis after conditioning the  CDC42  eQTLs and GWAS on the lead SNPs from an additional three independent  CDC42  eQTL signals. The  CDC42  association remained significant and sat just under the threshold for the HEIDI test ( P SMR  = 1.44 × 10 −6 ;  P HEIDI  = 1.11 × 10 −2 ).\nLike that performed for endometrial eQTLs, we also conducted an SMR-multi analysis using the eQTLGen data to check if heterogeneity results from multiple independent eQTL signals at this locus. Using all SNPs in each region resulted in an SMR-multi p-value ( P SMR-multi  = 7.45 × 10 −10 ) for  CDC42  similar to that observed from the standard SMR test ( Supplementary Fig. S7 ). However, we detected a second association between rs11801382 and long intergenic non-protein coding RNA 339 ( LINC00339 ) that failed to pass the standard SMR test, but passed the SMR-multi test ( P SMR-multi  = 2.73 × 10 −10 ,  P SMR  = 5.64 × 10 −2 ) ( Supplementary Fig. S7 ). The modified SMR test using the eQTLGen data found 18 associations that passed the SMR test ( P SMR  < 2.78 × 10 −4 ) including  LINC00339 ,  CDC42 ,  FGD6  and  VEZT  ( Table VIII ). However, none of the 18 pass the HEIDI test and therefore we could not distinguish between independent signals or one SNP affecting expression and disease risk in these regions.\nP-value significance thresholds:\nP SMR  < 2.78E-04\nP HEIDI  > 2.78E-03\nResults of modified SMR analysis using eQTLGen blood eQTLs and endometriosis GWAS summary statistics.\nEndometrial eQTLs also provide a valuable resource to functionally annotate genetic variants associated with other reproductive traits and diseases. We tested for overlap between independent Bonferroni significant eQTLs and SNPs associated with traits in the GWAS Catalog. We identified genetic variants that regulate endometrial gene expression and were associated with 288 traits and diseases ( P  < 5 × 10 −8 ) ( Supplementary Table SV ). This included various reproductive traits such as age of menarche onset, age of menopause onset, polycystic ovary syndrome, ovarian cancer and breast cancer. Two signals passed both the SMR and HEIDI test for age at menopause. These included loci on chromosomes 17 (rs2175957), associated with expression of Neighbor of BRCA1 LncRNA 2 ( NBR2 ) ( P SMR  = 2.27 × 10 −7 ,  P HEIDI  = 5.23 × 10 −2 ) and chromosome 20 (rs11699690), associated with expression of Copine 1 ( CPNE1 ) ( P SMR  = 8.16 × 10 −7 ,  P HEIDI  = 8.51 × 10 −2 ) ( Supplementary Fig. S8 ). Three signals passed both the SMR and HEIDI tests for epithelial ovarian cancer including two on chromosome 17; rs80028338 associated with expression of Leucine Rich Repeat Containing 37A ( LRRC37A ) ( P SMR  = 7.51 × 10 −10 ,  P HEIDI  = 3.44 × 10 −1 ) and rs17665188 associated with expression of Leucine Rich Repeat Containing 37 Member A2 ( LRRC37A2 ) ( P SMR  = 1.5 × 10 −9 ,  P HEIDI  = 6.68 × 10 −1 ) ( Supplementary Fig. S9a ). The third signal was located on chromosome 8; rs76837345 associated with expression of Charged Multivesicular Body Protein 4C ( CHMP4C ) ( P SMR  = 3.25 × 10 −6 ,  P HEIDI  = 9.73 × 10 −2 ) ( Supplementary Fig. S9b ). No signals passed the SMR test for epithelial endometrioid ovarian cancer.\n\nWe analysed genetic regulation of gene expression in the endometrium to determine how this relates to regulation in other human tissues, and whether genetic risk factors for endometriosis act through genetic effects on endometrial gene expression. The majority of common genetic effects on disease risk are located in non-coding regions of the genome and most likely act through regulation of gene expression in relevant pathogenic tissues ( Consortium  et al ., 2017 ;  Gamazon  et al ., 2018 ). We studied endometrial tissue because it is not represented in international projects like GTEx ( Consortium  et al ., 2017 ) and is a probable source of cells that initiate endometriosis lesions ( Sampson, 1927 ;  Anglesio  et al ., 2017 ;  Noë  et al ., 2018 ;  Suda  et al ., 2018 ). The transport of endometrial cells to the peritoneal cavity by retrograde menstruation as a cause of endometriosis was first proposed by  Sampson (1927)  and is supported by recent studies of somatic mutations in endometrium and endometriosis lesions ( Anglesio  et al ., 2017 ;  Noë  et al ., 2018 ;  Suda, et al., 2018 ). Our results show that the expression of genes located in genomic regions associated with endometriosis are significantly enriched in female reproductive tissues including uterus and endometrium, supporting this approach.\nWe first analysed genetic effects on gene expression (eQTLs) in endometrium and compared the eQTL profiles with tissues in GTEx ( Consortium  et al ., 2017 ) and large eQTL studies in blood. Generating an endometrial eQTL dataset from genotype and RNA-seq data, rather than our previous microarray data ( Fung  et al ., 2018 ) allowed a more accurate comparison with eQTLs reported in the GTEx ( Consortium  et al ., 2017 ). The large proportion (71.6%) of endometrial eQTLs also reported in the eQTLGen blood dataset highlights the potential power of using large datasets as a proxy for tissue shared eQTLs. Correlation in eQTL effects between endometrium and other tissues ranged from 0.54 to 0.72 with high correlations in eQTL effects between endometriosis and other reproductive tissues (vagina, uterus, breast, ovary). The relatively high correlations in eQTL effects is consistent with shared eQTL effects reported in the GTEx data ( Consortium  et al ., 2017 ;  Ongen  et al ., 2017 ), with more shared effects observed among tissues with greater biological similarity, for example among the reproductive tissues ( Consortium  et al ., 2017 ;  Ongen  et al ., 2017 ). Lower correlations between endometrium and blood in GTEx, endometrium and testis, and endometrium and brain is consistent with previous reports showing eQTLs in whole blood and testis have a high degree of tissue specificity ( Ongen  et al ., 2017 ).\nThe highest correlation of eQTL effects between endometrium and GTEx was observed with tissues of the digestive system. The underlying biology that leads to shared eQTLs between tissue is not yet clear. Digestive tract tissue is composed of an epithelial cell lining, endocrine epithelial glandular structures and mesenchymal derived support cells that function to secrete compounds required for tissue digestion and gut homeostasis ( Okumura and Takeda, 2017 ). The endometrium also has an epithelial lining with endocrine epithelial glandular structures and an endocrine secretory function designed to facilitate embryo implantation ( Hempstock  et al ., 2004 ), supported by mesenchymal derived stromal cells, surrounded by the smooth muscle of the myometrium. While many differences in the secreted substances exist, the molecular mechanisms to perform these roles may require similar gene regulation.\nIn contrast to the high correlation in cis-eQTL effects between tissues, no endometrial  trans -eGenes are reported in GTEx tissues and only three Bonferroni significant  trans -eGenes in endometrium are reported in the eQTLGen blood dataset ( Võsa  et al ., 2018 ). These findings are consistent with  trans -eQTLs being more tissue specific ( Grundberg  et al ., 2012 ;  Kirsten  et al ., 2015 ;  Consortium  et al ., 2017 ). The number of trans-genes identified in our study relative to sample size is similar to  trans -eQTL mapping for testis in GTEx which had the highest number of reported trans-eGenes ( n  = 35) ( Consortium  et al ., 2017 ). The larger number of trans-eQTLs in testis and endometrium may reflect the importance of  trans -SNPs in regulating gene expression in reproductive tissues. Larger eQTL studies in endometrium would be required to have sufficient power to accurately investigate  trans -acting genetic regulation of transcription.\nWe next analysed association between endometriosis and endometrial gene regulation. Differential expression analysis found no genome-wide significant differences in gene expression between endometriosis cases and controls following correction for multiple testing, consistent with our previous reports in eutopic endometrium ( Fung  et al ., 2017 ,  2018 ). An alternative approach combining gene expression and genotype data uses a powerful statistical method to impute gene expression from our eQTL data in a larger sample and conduct a TWAS. TWAS methods have been applied to identify functional loci in prostate cancer, obesity-related traits, Alzheimer’s, Crohn’s disease, diabetes and rheumatoid arthritis ( Gusev  et al ., 2016 ;  Mancuso  et al ., 2018 ;  Hu  et al ., 2019 ;  Nagpal  et al ., 2019 ). We conducted the first TWAS analysis for endometriosis using our endometrial gene expression data and identified 39 genomic regions associated with endometriosis. Five of the loci associated with endometriosis in the TWAS contain GWAS SNPs previously associated with the disease, including the  VEZT  locus,  KDR  locus,  GDAP1  locus and the  MLLT10  locus. The association between expression at the  VEZT  locus and endometriosis risk highlighted in the TWAS has been reported previously ( Sapkota  et al ., 2017 ;  Fung  et al ., 2018 ), and is supported by SMR analyses in this study. The other regions are novel. These have not been reported as genome-wide significant but all have nominal evidence of association in the GWAS studies. TWAS association does not imply causation but rather predicts differential gene expression between endometriosis cases and controls, with or without biological consequence. The correlation between predicted expression of genes at each genomic region limits the resolution to identify single target genes. Instead the method highlights potential candidates that may warrant further investigation.\nSeveral signals for genetic risk factors for endometriosis from GWAS ( Sapkota  et al ., 2017 ) and eQTLs map to the same regions of the genome. The overlapping signals are observed in different data sets and the overlap can occur by chance. Therefore, we used implementations of the SMR test to formally evaluate overlap in the signals and exclude chance overlap. We identified eQTLs in endometrium in which the same variant on chromosome 12 was associated with expression of both  VEZT  and  FGD6 , as well as risk of endometriosis. Several studies have reported associations between genetic variants at 12q22 locus near  VEZT  and increased risk of endometriosis ( Nyholt  et al ., 2012 ;  Rahmioglu  et al ., 2015 ;  Holdsworth-Carson  et al ., 2016 ;  Matalliotakis  et al ., 2017 ) and both genes were significantly associated with endometriosis in our TWAS analysis. The modified SMR analysis for expression in both endometrium and blood provides strong evidence that the same causal variant influences  VEZT  and  FGD6  expression and endometriosis risk. Both  VEZT  and  FGD6  play a role in plasma membrane, cell adhesion and cytoskeletal remodelling, all of which are important for development of endometriotic lesions ( Guo  et al ., 2011 ;  Holdsworth-Carson  et al ., 2016 ). Increased expression of  VEZT  has been reported in endometriosis cases, as has increased expression of epidermal growth factor receptor ( EFGR ) that is associated with expression of  FGD6  ( Ejskjær  et al ., 2009 ;  Meola  et al ., 2010 ;  de Graauw  et al ., 2014 ). Both genes have also been associated with expression of  CDC42 , another endometriosis risk gene ( Powell  et al ., 2016 ) responsible for cell division, growth and migration ( Miao  et al ., 2013 ;  Steenblock  et al ., 2014 ).\nWe observed a potential causal association between  CDC42  expression and endometriosis.  CDC42  passed the SMR test in blood and, following conditional analysis for multiple eQTL signals in the locus, came very close to passing the heterogeneity test. The multi-SNP SMR analysis identified a second variant associated with the expression of nearby  LINC00339  and endometriosis risk. These results are consistent with our previous studies where genetic regulation of  LINC00339  has been associated with endometriosis ( Powell  et al ., 2016 ;  Fung  et al ., 2018 ). Our results provide further support that genetic effects on endometriosis risk act through altered gene expression of  VEZT  and  FGD6  on chromosome 12 ( Holdsworth-Carson  et al ., 2016 ;  Powell  et al ., 2016 ) and through LINC00339 and possibly CDC42 on chromosome 1 ( Powell  et al ., 2016 ).\nUsing the standard SMR approach, we identified pleiotropic and potential causal associations between gene expression and age at menopause and ovarian cancer. The eQTLs for  NBR2  and  CPNE1  that were associated with age of menopause are reported in all 48 GTEx tissues but with varying ES.  NBR2,  which is located close to the tumour suppressor  BRAC1 , has been shown to have similar function to that of a tumour suppresser-regulating AMPK ( Liu  et al ., 2016 ). The largest genetic effect for  NBR2  was observed in ovary. We also identified eQTLs associated with ovarian cancer; one for  LRRC37 A was only reported in endometrium and the larger eQTLGen dataset with low expression across most GTEx tissues with the exception of testis. Interestingly, the eQTL for its paralog  LRRC37A2  was reported in 47 GTEx tissues, again with varying ES. Another signal passing the ovarian cancer SMR was for  CHMP4,  which is known as a prime candidate for epithelial ovarian cancer susceptibility due to its role in cell cycle regulation and regulation by TP53 ( Pharoah  et al ., 2013 ).\nOur study has important limitations. This is the largest study of genetic effects on gene expression in endometrium but it is small and lacks power in comparison with much larger eQTL studies in blood. Identifying and recruiting tissue donors is challenging due to the invasiveness of sampling, therefore limiting our ability to collect tissue on a large scale. Our analysis is conducted in fresh endometrial tissue consisting of multiple cell types. Consequently, expression levels are an average of expression from different cell types within the samples, which may mask smaller cell specific effects and identify only those large enough to be observed at a tissue level. Changes in cellular composition and cell activity across the cycle will contribute to variation in transcription across the cycle and between samples. We corrected for stage of the menstrual cycle in our analyses that will include changes in cellular composition across the cycle. Characterizing genetic regulation in individual cell types within the endometrium may also be important to understand the functional effects of disease risk variants. Future studies of expression in different cell types may identify novel cell-specific eQTLs if separation of the cells does not disrupt gene regulation and the studies have sufficient power. Techniques such as single-cell RNA-Seq offer an innovative solution to measure expression from individual cell populations; however, this technique would also introduce new practical, economic and computational challenges.\nIn conclusion, generation of an endometrial eQTL dataset using RNA-Seq and genome-wide genotyping data identified 327 novel genetic effects on transcription in endometrium. The ability to compare this dataset with publicly available eQTL datasets in GTEx and eQTLGen has identified high correlations in genetic effects between endometrium and both reproductive and digestive tissues and has allowed us to identify 68 endometrial  cis -eQTLs not observed in other tissues. Analysis of genetic effects on gene expression in endometrium provide further evidence that genetic risk factors for endometriosis act through expression of  VEZT ,  FGD6 ,  CDC42  and  LINC00339 . This was supported by the TWAS analysis with association for the  VEZT/FGD6  locus and endometriosis risk. The TWAS also identified a further 38 genomic regions harbouring potential target genes for functional follow-up. Expanding our knowledge of the genetic regulation in endometrium and integrating our data with publically available datasets creates an important resource to identify gene targets regulating female reproductive traits and diseases.\n\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.\nClick here for additional data file.","source_license":"CC0","license_restricted":false}