Integration of genome-wide association study and expression quantitative trait locus mapping for identification of endometriosis-associated genes

To determine whether genetic predisposition to endometriosis varies depending on ethnicity and in association with expression quantitative trait loci (eQTL) in a Taiwanese population. We conducted a genome-wide association study (GWAS) and replicated it in 259 individuals with laparoscopy-confirmed stage III or IV endometriosis (cases) and 171 women without endometriosis (controls). Their genomic DNA was extracted from blood and evaluated by the GWAS of Taiwan Biobank Array. Novel genetic variants that predispose individuals to endometriosis were identified using GWAS and replication, including rs10739199 (P = 6.75 × 10−5) and rs2025392 (P = 8.01 × 10−5) at chromosome 9, rs1998998 (P = 6.5 × 10−6) at chromosome 14, and rs6576560 (P = 9.7 × 10−6) at chromosome 15. After imputation, strong signals were exhibited by rs10822312 (P = 1.80 × 10−7) at chromosome 10, rs58991632 (P = 1.92 × 10−6) and rs2273422 (P = 2.42 × 10−6) at chromosome 20, and rs12566078 (P = 2.5 × 10−6) at chromosome 1. We used the Genotype-Tissue Expression (GTEx) database to observe eQTL. Among these SNPs, the cis-eQTL rs13126673 of inturned planar cell polarity protein (INTU) showed significant association with INTU expression (P = 5.1 × 10–33). Moreover, the eQTL analysis was performed on endometriotic tissues from women with endometriosis. The expression of INTU in 78 endometriotic tissue of women with endometriosis is associated with rs13126673 genotype (P = 0.034). To our knowledge, this is the first GWAS to link endometriosis and eQTL in a Taiwanese population.


Results
Assessment of population stratification. The demographic results are shown in Table 1. We performed a case-control GWAS to identify loci associated with increased risk of endometriosis in the Taiwanese population using an Affymetrix Axiom TWB array containing 653,291 SNP probes. We initially enrolled 126 endometriosis and 96 non-endometriosis controls from a Taiwanese population residing in Taiwan. After kinship analysis and strict quality control filtering, we analyzed 620,465 SNPs representing 95% of the array SNPs for the samples from the GWAS group. Multidimensional scaling analysis and results of permutation tests for identityby-state revealed no differences between the endometriosis and control groups, providing no evidence for strong population stratification (Fig. 1A,B). Quantile-quantile (Q-Q) plots were used to examine P value distributions, and the lambda value (λ) was 1.01 on the basis of the P value from the Cochrane-Armitage trend test, indicating no population substructure (Fig. 1C). In total, we found 33 top SNPs associated with endometriosis (P < 1 × 10 −4 ) (Supplementary Table 1).
GWAS and cross-platform validation. The GWAS analysis was conducted with 126 endometriosis cases and 96 controls. The Manhattan plot showed the result of genome-wide association analysis (−log 10 P) in chromosomal order for 620,465 SNPs test (Fig. 2). A minimum of 99% calling of Affymetric in both endometriosis cases and controls was selected for cross-platform validation using a Sequenom MassARRAY (Supplementary Table 1). Table 1) were tested in the replication stage and an independent cohort of 133 patients with endometriosis and 75 controls using a Sequenom MassARRAY and Q-PCR (Supplementary Table 2). After multiple test analyses using Bonferroni correction, the association was found to be not significant. In a combined analysis of the GWAS and replication cohorts, P values for 4 of the identified SNPs were lower than 10 −4 , which were not genome-wide significant (P < 5.0 × 10 -8 ) ( Table 2). We found that the SNPs rs10739199 (P = 6.75 × 10 −5 ) and rs2025392 (P = 8.01 × 10 -5 ) located at chromosome 9 in PTPRD (protein tyrosine phosphatase, receptor type D). Two SNPs (rs10739199 and rs2025392) were found in linkage disequilibrium (LD; D' = 0.961 and r 2 = 0.208, Fig. 3). After GWAS conditional analyses of these two SNPs, the P values were increased and indicated that they were only both associated. The P values for other 2 SNPs were lower than 10 -5 . These two SNPs were located at chromosome 14 (rs1998998, P = 6.5 × 10 −6 ) and at chromosome 15 (rs6576560, P = 9.7 × 10 −6 ). These were all replicated in the independent population and calculated in joint analysis ( Table 2). The P value of joint analysis is shown in Table 3 and Supplementary Table 3. However, the P value did not reach the standard genome-wide threshold (P value lower than 5 × 10 -8 ).
Genetic regulation of RNA transcription in endometrial tissue with defined genotypes-expression quantitative trait loci (eQTL) analysis. From discovery and replication studies of top 33 SNPs in trend tests, all the P values did not reach genome-wide significance. Expression quantitative trait loci (eQTL) explained the variation in expression levels of mRNA. To reveal the eQTL of the top 33 SNPs, we changed the significance threshold to 10 -3 in joint analysis and found 12 SNPs (Table 3). Among them, only rs13126673 is the putative cis-eQTL. The cis-eQTL analysis was performed to investigate potential association between the variants and expression levels of transcripts. The SNP rs13126673 was located at chromosome 4 in INTU. From Genotype-Tissue Expression (GTEx) project database v8, which contained the data of 322 testes samples from normal subjects, and it revealed that individuals with a CC genotype at rs13126673 have lower INTU expression compared with TT carriers, with a P value of 5.1 × 10 -33 (Fig. 7A). To further explore the eQTL in endometriotic tissues, 78 tissue samples from endometriosis patients with recorded SNP genotypes, were used for total RNA extraction and INTU expression was detected using RT-q-PCR. Of note, the C allele of SNP rs13126673 is the risk allele in our GWAS sample (OR =   ]. An eQTL analysis was performed to detect the effects of differing genotypes at SNP rs13126673 on the expression of the INTU transcripts; there was significant association between these genotypes and the expression of INTU transcripts observed (P = 0.034) (Fig. 7B).

Discussion
Several novel genetic variants of endometriosis were identified in this study, which, to our knowledge, represents the first report of a GWAS for endometriosis conducted in a Taiwanese population ( Figs. 1 and 2). In our twoindependent cohort, four novel loci for endometriosis were identified and replicated. Moreover, rs10739199 and rs2025392 were found in linkage disequilibrium (Fig. 3). After GWAS conditional analyses of these two SNPs, the P values were increased, and this indicated that they were only both associated. Imputation resulted in more strong signals (Figs. 4, 5, 6). rs13126673 was found that it the cis-eQTL of testes from the GTEx database and in endometriotic tissues (Fig. 7). The rs13126673 may affect the secondary RNA structure (Fig. 8). These results suggest a heritable component in endometriosis and provide new findings into the genetic risk factors of the disease. We found that the SNPs rs10739199 (P = 6.75 × 10 −5 ) and rs2025392 (P = 8.01 × 10 -5 ) located at chromosome 9 in PTPRD (protein tyrosine phosphatase, receptor type D) were associated with endometriosis. The other two SNPs were located at chromosome 14 (rs1998998, P = 6.5 × 10 −6 ) and at chromosome 15 (rs6576560, P = 9.7 × 10 −6 ). PTPRD is a member of the receptor protein tyrosine phosphatase (PTP) family. Mutations in PTPRD stimulate cell migration and growth in melanoma cell lines, enhance cell proliferation, and abrogate the dephosphorylation of Signal transducer and activator of transcription 3 (STAT3) in human astrocytes 15,16 . Several studies have shown that elevated STAT3 expression occurs in both endometriosis and endometrial www.nature.com/scientificreports/ cancer, suggesting STAT3 is a potential risk factor for both diseases [17][18][19][20] . In the present study, we found two endometriosis-associated SNPs, including rs10739199, and rs2025392 of chromosome 9 within the PTPRD gene. These novel genetic loci may provide new insights into the molecular basis of endometriosis. Moreover, deletions and mutations in PTPRD have been implicated in several tumor types, including endometrioid carcinomas in the Catalogue of Somatic Mutations in Cancer (COSMIC) database, and endometrial cancers 15 . Recently, a cross-disease GWAS meta-analysis revealed the rs2475335 SNP located within PTPRD to be associated with both endometriosis and endometrial cancer 21 . Further study of endometriosis and endometrial cancer models will be important to investigate the underlying biology of diseases-associated variants that increase the risk of both diseases 22 .
The major limitation of this study is the current sample size and it is too low to have real power to detect an association at a genome-wide level. However, endometriosis was reported to be common in patients with other benign gynecological diseases, especially uterine leiomyoma 23 . In the current study, all enrolled women had surgically confirmed diagnoses of endometriosis, and other benign gynecological diseases were used as the control group, hence the sample size was reduced. Recently, a GWAS of uterine myoma identified that eight novel genome-wide significant loci and four loci are also associated with endometriosis risk, suggesting overlapping genetic origins of uterine myoma and endometriosis 24 . In our study, we enrolled women with uterine myoma and other benign gynecological disorders as hospital-based controls rather than population-based controls, which does not inform the strength of endometriosis risk among other gynecological diseases. Because our controls were not healthy women, there may be some potential bias in our study. We also replicate the previous reported 14 SNPs and found that only rs13394619 was associated with endometriosis (P = 0.01968), this may due to different controls and populations. After multiple test analyses using Bonferroni correction, the association was found to be not significant. Further GWAS assessing potential susceptibility loci with genome-wide significance in different populations will verify genetic influences in the pathogenesis of endometriosis.
Moreover, we found that rs13126673 was at a novel locus on chromosome 4q28.1 with INTU and was associated with endometriosis. The rs13126673 SNP is an eQTL (expression quantitative trait loci) of INTU, which is a cilia and planar polarity effector with prominent ciliogenic roles in morphogenesis. rs13126673 is an intronic SNPs, which is not located in the predicted regulation regions of INTU transcription based on National Center for Biotechnology Information. Previous studies suggested that SNP may alter RNA/DNA structure and influence gene expression 25,26 . After prediction, the rs13126673 may change secondary RNA conformation (Fig. 8). The role of INTU in the pathogenesis of endometriosis will be worthwhile to study in the future.
In this study, we have provided the first genome-wide evidence in a Taiwanese population of four SNPs located in novel loci that were found to be associated with endometriosis. We have reported novel risk loci for the endometriosis-associated gene, PTPRD, that has been implicated in both endometriosis and endometrial cancer through cross-disease GWAS. www.nature.com/scientificreports/ Methods Study subjects. We recruited patients (n = 430) who underwent laparoscopic confirmation and histological analysis of endometriosis (compromising of 126 endometriosis cases from the GWAS and 96 cases from the discovery study) or who had no endometriosis (controls, comprising of 133 controls from the GWAS and 75 controls from the replication study) in the Center for Reproductive Medicine at Taipei Medical University Hospital. All diagnoses were confirmed by pathological analysis and divided into endometriosis and control groups. All women with endometriosis were diagnosed as stage III or IV endometriosis, the indication of control groups for laparoscopy includes myoma, fibrous adhesion, hydrosalpinx, ovarian cysts, teratoma, dermoid cysts, paratubal cysts, epithelial cysts, and simple cysts. A structured questionnaire was used by a trained researcher who interviewed subjects. Informed consent was obtained from each patient, and the study was approved by the Taipei Genotyping and quality control. Genomic DNA was extracted from blood using a DNA whole-blood kit, as per the manufacturer's instructions (Kurabo Industries, Osaka, Japan). The genotypes of each woman were performed by National Center for Genome Medicine (NCGM) at Academic Sinica using Axiom Genome-Wide TWB (Taiwan Biobank) Array Plate with a total of 653,291 SNPs. The Kinship analysis, quality control of genotypes for each SNP, total call rate, and Hardy-Weinberg Equilibrium (HWE) test were conducted by National     Table 1). SNP genotypes with a success rate over 97% and with over 99% concordance between the two platforms were then genotyped. The rs16911067 and rs16934324 were genotyping by Q-PCR using Taqman SNP genotyping assay and Taqman genotyping master mix (Thermo Fisher Scientific, MA, USA). The previous reported 14 SNPs were replicated using MALDI-TOP mass spectrometry. The multiple comparison correction was analyzed by Bonferroni test using the Graphpad Prism software (California, CA, USA).

Imputation.
To enhance coverage, the untyped SNPs were imputed by IMPUTE2 v2.3 using the 1000 Genome reference panel 27,28 . The National Center for Genome Medicine of Academia Sinica set up the haplotype inferences via the SHAPEIT method for optimizing the imputation rate 29 . We included a 500 kb buffer region on each side of the imputation region for elimination of edge effects and determined the uncertainty of imputed genotypes based on likelihood scoring in SNPTEST v2. Moreover, the frequentist association test of an additive model was used.
Expression quantitative trait loci (eQTL) analysis. The Genotype-Tissue Expression (GTEx) project database release the summary statistics of eQTL data of SNP rs13126673 for testis (n = 322) 30 . In our study, the eQTL was performed on the total of 78 tissue samples for endometriosis patients with recorded SNP genotypes. The expression of inturned planar cell polarity protein (INTU) was detected by reverse transcription quantitative polymerase chain reaction (RT-qPCR) and been described in detail elsewhere 31 . Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) mRNA was amplified using q-PCR with the forward primer 5′-gagtccactggcgtcttcac-3′ and reverse primer 5′-gttcacacccatgacgaaca-3′. The INTU mRNA was amplified using q-PCR with the forward primer 5′-tcagcgactcgggttcat-3′ and reverse primer 5′-cagccattcaggctcaaga-3′.
Prediction of RNA secondary structure. The rs13126673 SNP upstream and downstream DNA sequences (approximately 400 bp each) were retrieved by the dbSNP of the National Center for Biotechnology Information. The mfold program was used to predict the RNA structures of retrieved risk allele or normal allele of rs13126673 using default value and calculated the value of ΔG which represented the thermodynamic stability 12 . The smaller ΔG represented the more stable structure 12 .
Statistical analysis. The statistical method used for GWAS analysis is well-established by the National Center for Genome Medicine 32 . Possible population stratification could affect the association analysis and detection of this was performed using EIGENSTRA2.0. The variance inflation factor for genomic controls was also estimated. GC correction and genome-wide association were performed to compare allele and genotype frequencies between cases and controls using the Cochran-Armitage trend test. The P value distribution was showed in a quantile-quantile (Q-Q) plot. Adjustment for principle components suggested that inflation was not due to population stratification. The GWAS conditional analysis was performed as previously described 33 . The analysis of eQTL was performed with linear regression, using IBM SPSS statistics version 22 (New York, USA).