GWAS of five gynecologic diseases and cross-trait analysis in Japanese

We performed genome-wide association studies of five gynecologic diseases using data of 46,837 subjects (5236 uterine fibroid, 645 endometriosis, 647 ovarian cancer (OC), 909 uterine endometrial cancer (UEC), and 538 uterine cervical cancer (UCC) cases allowing overlaps, and 39,556 shared female controls) from Biobank Japan Project. We used the population-specific imputation reference panel (n = 3541), yielding 7,645,193 imputed variants. Analyses performed under logistic model, linear mixed model, and model incorporating correlations identified nine significant associations with three gynecologic diseases including four novel findings (rs79219469:C > T, LINC02183, P = 3.3 × 10−8 and rs567534295:C > T, BRCA1, P = 3.1 × 10−8 with OC, rs150806792:C > T, INS-IGF2, P = 4.9 × 10−8 and rs140991990:A > G, SOX9, P = 3.3 × 10−8 with UCC). Random-effect meta-analysis of the five GWASs correcting for the overlapping subjects suggested one novel shared risk locus (rs937380553:A > G, LOC730100, P = 2.0 × 10−8). Reverse regression analysis identified three additional novel associations (rs73494486:C > T, GABBR2, P = 4.8 × 10−8, rs145152209:A > G, SH3GL3/BNC1, P = 3.3 × 10−8, and rs147427629:G > A, LOC107985484, P = 3.8 × 10−8). Estimated heritability ranged from 0.026 for OC to 0.220 for endometriosis. Genetic correlations were relatively strong between OC and UEC, endometriosis and OC, and uterine fibroid and OC (rg > 0.79) compared with relatively weak correlations between UCC and the other four (rg = −0.08 ~ 0.25). We successfully identified genetic associations with gynecologic diseases in the Japanese population. Shared genetic effects among multiple related diseases may help understanding the pathophysiology.


Introduction
Uterine fibroma (UF), endometriosis, ovarian cancer (OC), uterine endometrial cancer (UEC), and uterine cervical cancer (UCC) are all common proliferative diseases arising from gynecologic organs. They are heterogeneous diseases with diverse range of proliferative and infiltrative properties. Clinical and epidemiological studies suggest that these diseases are mutually associated or often occur as comorbidity [1]. Studies of shared background risk, namely genetics, would offer understanding of the causes of these diseases, along with identifying targets to be treated.
For the past several years, genetic studies of gynecologic diseases have revealed only a limited number of significant associations [2][3][4][5][6][7][8][9][10]. Of note, the common risk genes wellknown from pedigree studies, such as BRCA1 and BRCA2 [11], have not been reported as either ovarian or breast cancer-susceptibility genes in the context of genome-wide association study (GWAS) [7,12]. This is largely because risk variants found in pedigree studies are usually rare among general population, which is unlikely detected in GWAS. So we focused on variants common among general population, including low frequency ones, to assess if these risk genes impose risks to general population, not just to specific families. The majority of the associations found in GWAS have small effect and polygenic nature. Detection of genetic associations in such cases depends on large sample size and ingenious analytical strategies. Genetic studies in the field of gynecology, including those reported from Biobank Japan Project (BBJ) [2,4,13], have been performed under the common logistic regression model. In this paper, to facilitate the detection of association signals and generalize the results, we conducted association studies under the liner mixed model (BOLT-LMM), with the largest GWAS data of Japanese population to date from BBJ. Advantages in adopting mixed models include that (i) they could account for both population stratification and cryptic relatedness, (ii) they could avoid confounding factors and provide robust association results, and (iii) they achieve increased statistical power for identifying genetic associations [14]. Large computational burden, a major limitation in mixed model methods especially when sample size is large, is overcome in BOLT-LMM [14]. In order to further enhance detection of novel loci, we also tried multiple association analytic approaches, where correlations between GWAS estimates among multiple related diseases are incorporated using MTAG [15]. Although MTAG utilizes summary level data and comparison between analyses based on row genotype data and summary level data might not be straightforward, the MTAG results were comparable to those from analyses under the common logistic model (mach2dat) and linear mixed model (BOLT-LMM) in a disease-specific manner. Since MTAG utilizes bivariate linkage disequilibrium (LD) score regression, where linear regression with liability threshold model is assumed and regression z-scores are assumed to follow standard normal distribution, which is different from the linear mixed model and non-normal distribution of regression z-scores assumed in BOLT-LMM, we utilized the results of mach2dat for applying MTAG. Also, we applied the reverse regression model using SCOPA [16], which utilizes raw genotype data, and produces estimates based on the best combination of phenotypes fitted to the model to obtain the maximized log-likelihood. We considered that SCOPA results could be comparable to those of joint analysis of all cases versus shared controls and random effect GWAS meta-analysis of different diseases in a multiple-disease-combined manner.
Clinically and epidemiologically, some of the gynecologic disorders harbor shared risk factors, such as age at menarche or menopause, and body mass index (BMI) [17][18][19]. These risk factors are at least partially influenced genetically and might be under the shared or pleiotropic effects of the genome [20]. To our knowledge, limited studies have investigated the shared genetic effects on gynecologic diseases, and they specifically look into relationship between two of the gynecologic diseases [21,22], relationship among histologic subtypes [23], or include gynecologic diseases as a part of multi-disease/trait study [24]. Only one of these studies analyzed the genetic correlation between specific pair of diseases under the linear mixed model [21]. To increase our understanding of shared genetic determinants influencing gynecologic diseases, here we report genetic correlations among the five gynecologic diseases in Japanese population using the linear mixed model approach.

Subjects
In total, 7315 cases with one of the clinically diagnosed five gynecologic diseases, including UC, endometriosis, OC, UEC, and UCC including cervical intraepithelial neoplasias, and 39,829 shared female controls without gynecologic diseases were enrolled from BBJ [13]. All the subjects provided written informed consent as approved by the ethical committee of RIKEN Yokohama Institute and the Institute of Medical Science, the University of Tokyo. This study was approved by the ethical committee of Osaka University Graduate School of Medicine. Females affected with nongynecologic malignancies and/or diseases thought to be strongly associated with the major histocompatibility complex (MHC) region were excluded from the control group. Related subjects were excluded in advance to avoid possible confounding. Carrier status of known risk genes such as BRCA1 and BRCA2, histopathological subtypes, and disease severity such as tumor sizes and clinical stages were not considered. Principal component analysis (PCA) was performed using EIGENSOFT (v6.1.4) Data manipulation was performed using PLINK software (v1.90b3.3).

Genotype imputation
Genotypes of the samples were obtained using either of the following genotyping arrays: (i) the Illumina Huma-nOmniExpressExome BeadChip or (ii) a combination of the Illumina HumanOmniExpress and HumanExome Bead-Chips. Genotype data is deposited on the Japanese Genotypephenotype Archive affiliated to the DNA Data Bank of Japan, via National Bioscience Database Center, Japan. The data are accessible on hum0014 at https://ddbj.nig.ac.jp/jga/ viewer/view/study/JGAS00000000114. For quality control (QC) of genotypes, we excluded variants meeting any of the following criteria: (i) call rate < 99%, (ii) P-value for departure from Hardy-Weinberg equilibrium (HWE) < 1.0 × 10 −6 , and (iii) number of heterozygotes less than five.
After we proceeded through these QC steps, we used Eagle (v2.3) for haplotype phasing without an external reference panel. We conducted whole-genome imputation using Minimac3 (v1.0.11) and the population-specific reference panel including multi-ethnic 2504 samples in 1000 Genomes Project (1KGP) phase 3v5a and deep whole genome sequencing of 1037 Japanese [25,26]. Variants with minor allele frequency (MAF) ≥ 1% in both case and control subjects, and imputation info r 2 ≥ 0.7 were selected for the following analyses.

Association analyses
Associations of the variants with each disease were separately evaluated (i) under the logistic regression model assuming additive effects of the allele dosages using mach2dat (v1.0.24), (ii) under the linear mixed model using BOLT(-LMM) (v2.2) [14], and (iii) by incorporating correlations between GWAS summary statistics using MTAG (v1.0.7) [15] with mach2dat results. Age, squared age, BMI, and the top 20 principal components (PCs) were included as covariates. PCs were calculated from 175,574 genotyped variants using EIGENSOFT (v6.1.4) under the same QC and pruning after excluding PCA outliers as described below in Heritability and genetic correlation section.
Joint analysis of all cases versus shared controls was executed under the logistic regression model using mach2dat, and was compared with random-effect metaanalysis of mach2dat results of each disease using RE2C (v1.04) [27] to correct for the overlapping samples and increase the power for detection. The same set of analyses was also conducted under the linear mixed model using BOLT-LMM and compared with the corresponding mach2dat results. We also applied the reverse regression model to the imputed genotypes and covariates-adjusted phenotypes using SCOPA (v1.0.14) [16]. The SCOPA results were compared with those from above analyses. Since RE2C* P-value, which is conditioned for the overlapping samples but is not conditioned for the fixed effect, was used for plotting, genome-wide significance threshold was defined as P = 5.0 × 10 −8 for all analyses [28]. Metasoft (v2.0.1) was applied to calculate metaanalysis heterogeneity index I 2 .

Candidate gene and functional annotation
For gynecologic risk loci outside of the MHC region, candidate genes were prioritized by FUMA [29]. We used the default settings except for the followings; we set the LD r 2 threshold to 0.5 to define LD structure of lead SNPs, reference panel population to 1KGP phase 3 EAS, and minimum allele frequency to 0.004. When there are no "mapped genes" then we listed the nearest coding or noncoding genes. Since genetic architecture within the MHC region is complex [30], it is displayed as "the MHC region" instead of the gene names. Functional annotations of the identified variants were also obtained by applying GAR-FIELD [31] and searching through HaploReg v4.1 (see URLs).

Heritability and genetic correlation
Heritability and genetic correlation were estimated using GCTA-GREML and phenotype-correlation-genotypecorrelation (PCGC-s) (v1.0.0) [32], respectively, with the same genotyped variant matrix, which passed the QC criteria as follows; we excluded samples with call rate < 99%, variants with call rate < 99%, variants with MAF < 5%, variants with HWE P-value < 1.0 × 10 −6 , variants in the MHC with flanking region spanning from 24 to 34 Mb of chromosome 6, then pruned variants at LD r 2 cutoff of 0.50. Disease prevalence was defined as described elsewhere [1,33,34]. Age, squared age, BMI, and 20 PCs were incorporated as covariates. BOLT-REML failed to converge in analyzing heritability, and GCTA-GREML failed to converge in analyzing genetic correlation, therefore we adopted PCGC-s that would definitely yield the results through Haseman-Elston regression model. PCGC-s does not produce standard errors or p-values for the estimates.

Subjects
Numbers of the subjects eligible for each GWAS were as follows; 5236 for UF, 645 for endometriosis, 647 for OC, 909 for UEC, and 538 for UCC cases where those who have multiple diseases were allowed to enroll in each corresponding GWAS, and 39,556 shared female controls.
Characteristics of the samples included in the GWASs of the five gynecologic diseases are shown in Table 1.

Genotype imputation
Whole-genome imputation yielded 61,608,817 variants in total. Of these variants, 7,645,193 passed the quality control as described above (MAF ≧ 1% in both cases and controls, and imputation info r 2 ≧ 0.70).
There were five associations that reached genome-wide significance by BOLT-LMM but not by mach2dat; rs12415148:T > C at 10q24 for UF, rs79219469:C > T at 16q12 and rs567534295:C > T at 17q21.31 for OC, and rs150806792:C > T at 11p15 and rs140991990:A > G at 17q24 for UCC. Of these, four were low-frequency variants (1% < MAF < 5%). Although some associations were also detected as significant by MTAG, there were no variants that surpassed the genome-wide significance only by MTAG.
Joint analysis of all the cases and controls, randomeffect meta-analysis of single disease GWAS, and reverse regression analysis We first performed two additional analyses; (i) joint analysis of all the combined cases of the five gynecologic diseases versus shared controls and (ii) random-effect metaanalysis of five single-disease association studies correcting for the overlapping samples using RE2C.  Joint analysis of all cases versus shared controls yielded four associations that surpassed the significance threshold (rs7412010:G > C at 1p36, rs12225799:C > G at 11p15, rs17332320:G > T at 22q13, and the MHC region at 6p22, Fig. 2 and Table 3). Out of the four associations, three (rs7412010:G > C at 1p36, rs12225799:C > G at 11p15, and rs17332320:G > T at 22q13) were identical to those detected in single disease analysis of UF, and were detected both by mach2dat and BOLT-LMM, while the remaining one, the MHC region, was identical to that found in single disease analysis of UCC, and was significant only in the analysis by BOLT-LMM. For all associations that surpassed the genome-wide significance, BOLT-LMM showed more significance than mach2dat, as was shown in the singledisease associations.
We next performed reverse regression analysis using SCOPA. When the MHC region is counted as one locus, we BOLT-LMM mach2dat   The underlined odds ratios (OR) and P-values indicate the analytic method which showed the most significant association for each SNP among the four methods, joint analysis of all cases versus shared controls and random-effect meta-analysis using mach2dat and BOLT-LMM, respectively OR in BOLT-LMM is adjusted using case fraction u, using the formula; log(OR) identified 12 loci in association with gynecologic diseases (Fig. 2 and Table 4). Of these, seven loci were best explained when the multiple gynecologic diseases were combined in the regression model. In addition to the three novel loci identified in single disease GWAS and meta-analysis, three loci (rs73494486, chr9:g.101341851 C > T at 9q22, GABBR2, P = 4.8 × 10 −8 ; rs145152209, chr15:g.84077212 A > G at 15q25, SH3GL3/BNC1, P = 3.3 × 10 −8 ; and rs147427629, chr21:g.40419321 G > A at 21q22, LOC107985484, P = 3.8 × 10 −8 ) were novel findings. The top SNPs in the previously identified loci were almost the same as those detected in single disease GWAS and meta-analysis (Tables 2-4). Annotations obtained by searching through HaploReg v4.1 are summarized in Supplementary Table 12.
When comparing the five patterns of the analyses, namely, combination of two conjoining methods (joint analysis of all cases versus shared controls, and randomeffect meta-analysis of the five GWAS correcting for the overlapping subjects), and two analytic methods (the usual logistic model using mach2dat, and the linear mixed model using BOLT-LMM), and reverse regression model, the association results were heterogeneous. These results would suggest that applying multiple methods would offer more opportunity to detect significantly associated loci.

MHC region
Within those detected loci, the variants showing the most significant associations are almost identical among the analytic methods. However, the MHC region at chromosome 6, which surpassed genome-wide significance in single-disease GWAS of UCC, joint analysis, and randomeffect meta-analysis, showed the most significant associations at different variants among the methods as previously suggested [30]. Fine-mapping and identification of causal variants of the MHC region by high resolution HLA imputation would be warranted [30,35].

BRCA1
BRCA1 is a tumor suppressor gene well-known from the pedigree studies of familial breast and ovarian cancers. In this study, we identified the significant association with OC at this gene locus in the GWAS approach for the first time. While previous genetic studies mostly focused on coding variants of BRCA1 [36], the associated variant, rs567534295:C > T, was the low-frequency noncoding variant (1% < MAF among controls < 5%), located within the intron between exon 22 and exon 23, and all the variants Based on hg19 c Reference (Ref) and alternative (Alt) alleles on forward strand d Alternative allele frequency among control subjects e SCOPA produces P-value of the best model in moderate LD (r 2 > 0.5) with rs567534295:C > T were not located within the coding region. Rs567534295:C > T is monomorphic in non-east Asian populations (in 1KGP phase3v5a), suggesting population-specific risk of the variant on OC. To make functional annotation of this lowfrequency risk variant, we assessed the existing expression quantitative trait locus (eQTL) databases (GTEx, HGVD, Ishigaki et al. [37]) and applied artificial intelligence(AI)based prediction algorithm on variant functions (ExPecto [38]). However, known eQTL databases and annotation tools do not contain these variants, and AI-based prediction showed that BRCA1 expression was not affected by the variant rs567534295:C > T, suggesting the contribution of undetermined mechanisms to the pathophysiology of OC.

Heritability and genetic correlation
The disease prevalence used for estimating heritability is shown in Table 1. GCTA-GREML applied to the genotype matrix of each disease produced very stringent results; heritability estimated from genotyped SNPs (h 2 SNP in Table 1) ranged from 0.0260 for ovarian cancer to 0.220 for endometriosis (Table 1). We observed that the h 2 SNP is about the fifteenth to the two thirds of the disease heritability reported in epidemiological studies (h 2 g in Table 1). Genetic correlations among the five gynecologic diseases under the linear mixed model, calculated directly from the genotyped SNPs, are shown in Fig. 3. All correlations among the four gynecologic diseases except for UCC were directionally positive, and stronger correlations were detected between endometriosis and OC (r g ≥ 1.00), OC and UEC (r g ≥ 1.00), and UF and OC (r g = 0.79). While relatively weaker and sometimes negative correlations were depicted between UCC and the four other gynecologic diseases (r g = −0.08-0.25, Fig. 3).

Discussion
In this study, we detected nine significant associations with three of the five gynecologic diseases. Four out of the nine associated loci were identified by applying the linear mixed model but not by the usual logistic regression model. For example, rs12415148:T > C at STN1/SLK associated with UF is a known risk locus for UF but was not detected by the usual logistic regression model approach with the current sample size. This supports the application of the linear mixed model approach to the case-control studies to enhance the detection of genetic associations.
We identified two novel loci associated with OC, rs79219469:C > T and rs567534295:C > T. We compared our results with the publicly available summary statistics from the largest European OC GWAS. We found that chr16 variant rs79219469:C > T exists in the summary but its effect allele (T) frequency is only 0.0046, with its association p-value 0.84. We also looked for the chr17 variant rs567534295:C > T, however, this does not exist in the available summary (monomorphic in 1KG Europeans). We considered that disease risk of these variants were rather population-specific in Japanese (and east Asians).
We identified two novel loci associated with UCC, rs150806792:C > T and rs140991990:A > G. Rs150806792: C > T at INS-IGF2 locus is suggested to be associated in the pathogenesis of UCC through the activation of insulin-like growth factor pathway, as reported in the cases of colorectal tumors [39]. Rs140991990:A > G is located at SOX9 locus, which is a member of SOX family and its family member SOX14 is reported to be involved in p53 signaling pathway in a UCC cell line [40]. Also, inhibition of SOX9 is reported to increase radiosensitivity in gastrointestinal cancer [41]. These suggest the roll of SOX9 in the pathogenesis of UCC.
Of the nine identified associations, three top variants; rs567534295:C > T for OC, and rs150806792:C > T and rs140991990:A > G for UCC, are variants found only in Japanese and/or East Asians with rare to low-frequency but not observed in other populations, according to 1KGP phase 3v5a data. These three were also the variants better detected under the mixed model. Of note, the well-known BRCA1 locus, encompassing the low frequency noncoding variant of rs567534295:C > T, was shown to be significantly associated with OC for the first time in the context of  Fig. 3 Cross-trait evaluation of genetic correlation among five gynecologic diseases. Genetic correlations among five gynecologic diseases calculated under the linear mixed model by Haseman-Elston regression using PCGC-s. Correlation is expressed by the color and size of square on the right upper triangle, while represented in digits on the left lower triangle. Asterisks indicate that the real output value exceeded one but was set to one for display purpose GWAS [7]. We note that risk identification of rs567534295: C > T could have been achieved only by using the large population-specific imputation reference panel [25,26], with adequate GWAS sample size, and by using ingenious analytical methods. While the GARFIELD analysis suggested functional annotations of genetics of the phenotypes, cell specificity observed in the results were still relatively nonspecific. Further approaches would be warranted to further elucidate undetermined disease mechanisms.
In the joint analysis, since sample sizes are different among the diseases, signals are most likely to be driven by the disease with the largest samples. As expected, three signals were concordant with those detected in UF GWAS, having the largest samples among the five GWASs, however, signals in the MHC region were also significant under the linear mixed model, which were concordant with those detected in UCC GWAS, having the smallest samples. This showed that strength of the signals in joint analysis can be either amplified or attenuated depending on the sample sizes and correlations among the diseases at the signal of the interest.
A novel association was also discovered by random-effect meta-analysis of the five GWASs adjusting for the overlapping samples. The top associated SNP rs937380553:A > G at 2p16 locus is located within noncoding LOC730100 gene. The function of this lncRNA is not well investigated and further studies are required to elucidate its contribution to the pathogenesis of gynecologic diseases, especially, possibly shared effect among endometriosis, OC, and UEC. We also applied a novel approach, MTAG, which incorporates correlations among multiple GWAS estimates to enhance detection, to gynecologic diseases. In some, but not all, MTAG results were the most significant among those of the usual logistic regression model and the linear mixed model. Therefore, by applying multiple analytic methods, we can have more opportunity to identify novel associations.
In the reverse regression analysis, we identified three additional novel associations. Rs73494486:C > T at GABBR2 locus is associated with combined phenotypes of UF and OC. This locus has neither been reported in the GWAS of UF or OC, however, GABBR2 is suggested to have an important role in EGFR signaling through the ERK1/2 pathway, as reported in lung adenocarcinoma [42]. Rs145152209:A > G at SH3GL3/BNC1 locus is associated with UF. SH3GL3 and BNC1 are both neither reported in the context of UF, however, SH3GL3 is reported as a colorectal cancer-associated gene [43], and BNC1 is reported to have association with pancreatic cancer [44,45]. These suggest that this locus may play an role in proliferative property of UF in some undetermined mechanisms. Rs147427629:G > A at LOC107985484 locus is associated with OC + UEC. This noncoding RNA is not studied well, however, the similarity of OC and UEC supports the mutual pathophysiology conferred by this non-coding RNA.
We showed the genetic correlations under the mixed model among the five gynecologic diseases. As expected from the past genetic and epidemiologic studies [21][22][23], these correlations were mostly directionally positive. The correlation was stronger between OC and UEC, and between endometriosis and OC, compared with those between UCC and the other four gynecologic diseases, which was concordant with the epidemiological findings [1,[20][21][22]. Thinking of the pathophysiology, where UCC is mainly caused by infection of human papilloma viruses, which is distinct from other gynecologic diseases, and where some histological subtypes of OC very often co-occur with endometriosis [1], this result would be considered rational. Although genome-wide additive effects are correlated given the strong genetic correlation and similarity of SNP effect sizes among these gynecologic diseases, we could detect only one additional shared locus in the cross-trait meta-analysis. This was because many diseases have polygenic nature and most of the variants have too small effect sizes to be detected. By increasing sample sizes, more shared loci with small effects are expected to be detected. This encourages cross-trait meta-analysis to support improved power to detect shared loci as described above.
Strength of our study includes the followings: First, we conducted the GWASs with the largest sample sizes among Japanese population, which facilitated the detection of novel associations. Second, we performed association analyses under the linear mixed model and also applied MTAG. We demonstrated that application of the mixed model and/or integrating correlations among multiple diseases increased the power of detection. Third, we used the population-specific imputation reference panel, which contributed to the increased number of variants with higher accuracy, especially those specific to the Japanese population such as the risk noncoding variant at BRCA1. Fourth, we demonstrated that combining the GWASs of multiple diseases in a random-effect metaanalysis revealed a novel candidate association. Fifth, we revealed genetic correlations among multiple gynecologic diseases. Shared genetic etiology encourages the investigation of common pathophysiology of the related diseases.
Although our study is the largest GWAS in the Japanese population, our limitation includes the lack of replication study. This is the task for our future study.
In conclusion, we successfully identified nine significant genetic associations with three gynecologic diseases including four novel ones, by applying association analysis under the mixed model and incorporating correlations among multiple GWAS estimates. Further, cross-trait meta-analysis identified five loci including one novel association which is suggested to be a shared risk locus. We also disclosed genetic correlations among multiple gynecologic diseases. We propose to apply new methodologies to increase detection power, and crosstrait analysis to assess shared risks. ment (AMED). This study was supported by Bioinformatics Initiative Author contributions TM, YU, KT, YK, and YO designed the study. TM, SKL, and MA performed statistical analyses. MH, KM, YM, and MK collected and managed the samples. MK conducted genotyping. TM and YO wrote the manuscript. All authors helped interpret the data, reviewed, and edited the final manuscript, and approved the submission.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval This study was approved by the ethical committee of Osaka University Graduate School of Medicine.
Informed consent All the subjects provided written informed consent as approved by the ethical committee of RIKEN Yokohama Institute and the Institute of Medical Science, the University of Tokyo.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.