A common variant within the HNF1B gene is associated with overall survival of multiple myeloma patients: Results from the IMMEnSE consortium and meta-analysis

Diabetogenic single nucleotide polymorphisms (SNPs) have recently been associated with multiple myeloma (MM) risk but their impact on overall survival (OS) of MM patients has not been analysed yet. In order to investigate the impact of 58 GWAS-identified variants for type 2 diabetes (T2D) on OS of patients with MM, we analysed genotyping data of 936 MM patients collected by the International Multiple Myeloma rESEarch (IMMENSE) consortium and an independent set of 700 MM patients recruited by the University Clinic of Heidelberg. A meta-analysis of the cox regression results of the two sets showed that rs7501939 located in the HNF1B gene negatively impacted OS (HRRec= 1.44, 95% CI = 1.18–1.76, P = 0.0001). The meta-analysis also showed a noteworthy gender-specific association of the SLC30A8rs13266634 SNP with OS. The presence of each additional copy of the minor allele at rs13266634 was associated with poor OS in men whereas no association was seen in women (HRMen-Add = 1.32, 95% CI 1.13–1.54, P = 0.0003). In conclusion, these data suggest that the HNF1Brs7501939 SNP confers poor OS in patients with MM and that a SNP in SLC30A8 affect OS in men.


INTRODUCTION
Multiple myeloma is an incurable and heterogeneous plasma cell neoplasm that affects about 6.3 per 100.000 people per year worldwide (i.e., 25.850 new cases only in 2015) and represents 1.6% of all cancers and 2% of all cancer deaths [1]. In spite of the widespread use of proteasome inhibitors and immunomodulatory drugs, which has dramatically improved the life expectancy of MM patients over the last few decades [2,3], MM survival still remains poor with a 5-year survival of 46.6% (SEER Cancer Statistics Review, http://seer.cancer.gov/ csr/1975_2012/).
Epidemiological and observational studies have consistently identified several factors that affect MM patient survival such as age at diagnosis [4,5], stage at diagnosis (coded by either the Durie-Salmon staging system (DSS) [6] or the International Staging System (ISS)) [7], Eastern Cooperative Oncology Group (ECOG) performance status [8], renal failure [9,10], high plasma cell proliferative rate [11,12], high lactate deshydrogenase (LDH) levels [13] and chromosomal abnormalities [14][15][16][17][18]. Increasing evidences point towards a positive correlation of pre-existing type 2 diabetes (T2D) with MM risk [19] but also with the appearance of severe clinical complications [20][21][22][23] and patient survival [24,25]. In this regard, Chiu et al. (2006) reported that high level of postload glucose was associated with increased risk of mortality in hematological malignancies [24] whereas Chou et al. (2012) reported that MM patients with preexisting T2D have 50% higher all-cause mortality compared with non-diabetic patients [25]. These observations might be explained, at least in part, by the stimulatory effects of T2D-associated hyperglycaemia, insulin resistance and resulting hyperinsulinemia on MM cell growth [26,27] but also by the deregulation of tumour-suppressor genes linked to T2D (such as CDKN2A-2B, KCNQ1, HNF1B) [28][29][30] that might lead to uncontrolled cell proliferation, cell differentiation and disease progression and, consequently, to shorter survival periods. In support of this notion, CDKN2A-2B genes have been found to be frequently hypermethylated in MM [31][32][33] whereas loss of expression of KCNQ1 has been associated with poor overall survival in cancer patients [34]. Furthermore, emerging evidences also suggest that the activation of certain T2D-related genes (such as NOTCH2) may induce MM cell migration from the infiltrated site to different bone marrow districts [35] and promote osteoclast formation [36], which is a process intimately related to proliferation and long-term survival of MM cells [37].
Although germline variants may influence the susceptibility of MM [38][39][40][41][42][43] and survival [39,[44][45][46], the knowledge regarding the role of diabetogenic variants in modulating the risk of MM and survival remains scarce. We have recently reported that diabetogenic variants influence MM risk [47] and recent genome-wide association studies (GWAS) have also suggested the involvement of genetic variants within the MTHFD1L, AKAP12 and FOPNL loci in determining MM patient survival [39,44,45] but also an indirect implication of diabetogenic genes such as TCF7L2 [45]. Johnson et al. (2016) reported in their GWAS a strong association of rs12374648, which maps to a binding site for the transcription factor TCF7L2, with MM overall survival (OS) and proposed a functional mechanism of this variant to modulate the synthesis of purines and the regulation of cell cycle [45].
Based on these findings, we explored for the first time the relationship between diabetogenic variants and OS of MM patients in a study developed in the context of the International Multiple Myeloma rESEarch (IMMENSE) consortium. We attempted to confirm our findings by analysing GWAS data on an independent set of German MM patients (Heidelberg cohort) [45].

RESULTS
The demographic and clinical characteristics of the MM patients included in the IMMENSE (n = 939) and Heidelberg (n = 700) cohorts are listed in Table 1 All SNPs tested showed genotype frequencies consistent with the HWE (P > 0.001) and the observed allele frequency for all selected SNPs was in accordance with Hapmap data. When we evaluated the effect of selected polymorphisms on MM OS in the IMMENSE population, we found that 6 SNPs showed a noteworthy association with OS. The most relevant effect was observed for the HNF1B rs7501939 SNP that was associated with poor OS when recessive and log-additive models of inheritance were assumed (HR REC = 1.49, 95% CI 1.11-2.00, P = 0.008 and HR ADD = 1.34 95% CI 1.13-1.59, P = 0.001, respectively; Table 2). Patients harbouring the nondiabetogenic HNF1B rs7501939T/T genotype showed a median survival time (MST) significantly shorter than those carrying the C allele (MST T/T = 81.91 months vs. MST C/C+C/T = 101.42 months; Figure 1A). This result was confirmed with the 700 MM patients recruited from the University Clinic of Heidelberg (HR REC = 1.40, 95% CI 1.06-1.84 and MST T/T = 74.4 months vs. MST C/C+C/T = 97.2 months; Table 3 and Figure 1B). The result of the metaanalysis for this SNP remained significant after correction for multiple testing (HR Meta-Rec = 1.44, 95% CI 1.18-1.76, P = 0.0001, I 2 = 0.0% , P Het = 0.74; Table 3). According to publicly available eQTL data for human peripheral blood mononuclear cells, the risk allele (T) was associated with higher HNF1B mRNA expression levels (Z score = 3.31, P = 9.23•10 -4 and FDR = 0.23). However, we could not validate this finding using eQTL data on plasma cells from 665 German MM patients (P = 0.60; Supplementary  Table S1). Nonetheless, according to Haploreg and ENCODE annotation data, the rs7501939 SNP resides near of a poised promoter in many cell lines including a lymphoblastoid and human stem cell lines (GM12878 and H1-HSCs) that might be rapidly activated upon specific stimuli. In addition, this SNP was predicted to change binding motifs for 2 regulatory transcription factors (CEBPB and p300) and mapped among enhancer histone marks in primary naïve and memory forms of cytotoxic T cells (CD 8 + ) and helper T cells (CD 4 + ) from peripheral blood.
We also observed significant associations at P < 0.05 for SNPs within CDKN2A-2B, GCKR, KCNQ1 and SLC30A8 genes with OS in the IMMENSE population. Thus, patients carrying the KCNQ1 rs2074196T and SLC30A8 rs13266634T alleles or the GCKR rs1260326T/T genotype had an increased risk of death whereas subjects bearing the CDKN2A-2B rs564398C/C genotype showed longer OS (Table 2 and Supplementary Tables S2-S4). The association of the SLC30A8 rs13266634T allele with OS was confirmed in the Heidelberg population and the metaanalysis showed that the presence of each additional copy of the SLC30A8 rs13266634T allele was associated with poor OS (HR Meta-Add = 1.22, 95% CI 1.09-1.37; Table 3). Although the association of the SLC2A2 rs11920090 SNP with OS was not significant in the IMMENSE population, we observed a significant association of this variant with MM survival in the Heidelberg population that remained significant in the pooled analysis. Patients harbouring the SLC2A2 rs11920090T allele showed a better survival compared with those carrying the A/A genotype (HR Meta-Dom = 0.80, 95% CI 0.66-0.97; Table 3). The meta-analysis also showed a weak association of the TCF7L2 rs7903146T allele with better survival that was neither significant in the IMMENSE population nor in the Heidelberg cohort (HR Meta-Dom = 0.84, 95% CI 0.71-0.98). Based on Haploreg data, the missense rs13266634 SNP was predicted to change binding motifs for transcription factors implicated in tumorigenesis (AP1 and PAX5) and mapped on enhancer histone marks in several human embrionic stem cell lines. In addition, this polymorphism affects binding to 5 proteins implicated in cancer development (CCNT2, GATA2, TAL1, KAP1 and CTCF). On the other hand, the rs11920090 and rs7903146 SNPs mapped among enhancer and promoter histone marks in bone marrow-and/or adipose-derived mesenchymal stem cells. In addition, the rs7903146 SNP was predicted to alter the binding site of 7 transcription factors. However, despite the consistency and the potential interest of these findings, none of the associations of the SLC2A2, SLC30A8, and TCF7L2 SNPs with OS remained significant after correction for multiple testing and, therefore, require further confirmation. Given the lack of genetic information in the GWAS conducted in the Heidelberg population for SNPs within MADD and KCNQ1 genes, we imputed genotypes to test whether the preliminary associations observed in the IMMENSE population could be validated. Although there was no imputed data available for the KCNQ1 SNP, the metaanalysis of IMMENSE data with imputed genotypes of MADD variants in the Heidelberg cohort suggested a link between this locus and MM survival (HR Meta-Rec = 0.75, 95% CI 0.57-0.99; Table 3).
Based on the evidences that point toward the existence of gender-associated differences in survival for patients with MM [48], we decided to carry out a gender-stratified analysis. This analysis revealed gender-specific associations for SNPs within or near the ADAMTS9, KCNJ11, PROX1 and SLC30A8 genes with OS. We found that men carrying the KCNJ11 rs5215C or SLC30A8 rs13266634T alleles had poorer OS compared with those harbouring the wild type genotype whereas an opposite but not significant effect was seen in women (P Interaction = 0.022 and P Interaction = 0.057, respectively; Table 2 and Supplementary Table S2-S4). We also observed that women carrying the PROX1 rs340874G allele or the ADAMTS9 rs4607103T/T genotype experienced a poorer survival with an opposite but not significant effect in men (P Interaction = 0.016 and P Interaction = 0.024). In order to confirm these gender-specific associations, we performed a   Table S5). This gender-specific association remained significant at the experiment-wide significance threshold. On the other hand, although it was not statistically significant in the analysis of the IMMENSE population, the pooled analysis also showed that men carrying the BCL11A rs10490072C allele had a poorer OS compared with those carrying the wild type genotype whereas no effect was seen in women (HR Men = 1.37, 95% CI 1.10-1.70). Finally, we observed in the pooled analysis that women bearing the PRC1 rs8042680A allele or men carrying the PROX1 rs340874G allele or the MADD rs7944584T/T rs1531343 RPSAP52, HMGA2   ) and regulatory T cells from peripheral blood.

Previous
population-based studies have demonstrated the impact of GWAS-identified variants for T2D on cancer susceptibility [47,[49][50][51][52] and patient survival [53]. However, despite these important research advances, there is still a noticeable lack of information regarding the role of T2D-related variants in modulating patient survival especially in hematological malignancies. In this scenario, we decided to investigate for the first time to our knowledge the relationship between 58 genetic variants associated with T2D identified by GWAS and OS of MM patients.
The analysis of the IMMENSE consortium data revealed a significant association of the intronic HNF1B rs7501939 SNP with poor OS. We successfully replicated this association in a large and independent population recruited by the University Clinic of Heidelberg. However, although a positive correlation between this variant and eQTL data on PBMCs has been reported [54], we failed to find correlation between the risk allele and HNF1B mRNA expression levels on plasma cells from a large cohort of MM patients. This suggested that the effect of this variant on overall survival is not mediated by changes in transcriptional activity of the  gene. Nonetheless, given that HFN1B contains multiple independent SNPs or haplotypes that have been associated with HNF1B mRNA expression [55,56] and methylation [57] levels but also with the risk of developing several types of cancer [55][56][57], it seems to be reasonable to consider the possibility that other SNPs within this locus and showing a stronger association with OS could explain better the link between the HNF1B and clinical outcome. However, when we analysed imputed common SNPs from the GWAS conducted in the Heidelberg cohort, we could not find any stronger association signals with OS in the region, which suggested that the HNF1B rs7501939 SNP or perhaps a rare SNP in LD with it might be responsible of the observed effect. Future fine-mapping studies encompassing common but also rare variants within or near the HNF1B locus are needed to elucidate whether a rare variant or haplotype might account for the observed effect.
HNF1B contains 9 exons and expands over 58 kb on chromosome 17p21 [58,59]. It encodes for a transcription factor that has been associated with multiple clinical features including early-onset of T2D [56]. In line with this, it has been also suggested that HNF1B may induce impaired glucose tolerance and attenuated insulin sensitivity in a miRNA-dependent manner [60], which might lead to an enhanced insulin secretion and the activation of the IGF1 pathway, an important factor mediating myeloma cell growth, proliferation and cell maturation [27,59,61]. Alternatively, it has been postulated that HNF1B is able to influence cancer cell survival by promoting the activation of NFkB pathway or through the inhibition of mitochondria-associated apoptotic signals [62]. In support of the tumorogenic effect of HFN1B, it has also been reported that it may act as an oncogene [63] and that the HNF1B gene is amplified in 23% of all cancers and in about 5% of all haematological malignancies (http://broadinstitute.org/ tumorscape). On the contrary, it has also been reported that HFN1B may act as tumour suppressor gene [56] and that its expression may largely vary depending on the target tissue. Whereas HNF1B has been found to be overexpressed in ovarian clear cell carcinomas [58] and prostate [64] or endometrial [65] cancers and its silencing induces apoptosis of cancer cells [58], it has been found to be down-regulated in serous epithelial ovarian cancer [57] and colorectal, gastric and pancreatic cancer cell lines [66]. In addition, it has been reported that the downregulation of HFN1B gene is associated with progression in hepatocarcinoma [67] and poor prognosis in renal [68] and prostate [69] cancers. Considering all the above but also the fact that the association of the HNF1Brs7501939 SNP with OS was driven by a non-diabetogenic (T) allele that does not affect HNF1B mRNA expression, we hypothesize that the effect of this variant to contribute to tumour progression in MM might be mediated by a noninsulin-dependent mechanism. There was a reasonable amount of regulatory data for the HNF1Brs7501939 SNP that supported evidence of the active role of the HNF1B locus. However, whether elevated HNF1B levels lead to tumour transformation and disease progression is not yet understood and functional studies to examine whether HNF1B variants influence cancer prognosis are lacking.
Another interesting finding of this study was the association of the BCL11A, MADD, PRC1, PROX1, SCL30A8, SLC2A2 and TCF7L2 SNPs with OS. We found an overall association of the SLC2A2 rs11920090T and TCF7L2 rs7903146T alleles with better OS whereas the association of the SLC30A8 rs13266634T , BCL11A rs10490072C , PRC1 rs8042680A and PROX1 rs340874G alleles or the MADD rs7944584T/T genotype with OS was restricted to male or female genders. Despite the potential interest of the associations observed for these SNPs with OS, only the association of the SLC30A8 rs13266634 SNP with poor OS in men reached significance at experiment-wide significance threshold. This result suggested a key role of the SLC30A8 locus in the modulation of overall survival. However, given the consistency of the overall or genderspecific associations observed for BCL11A, MADD, PRC1, PROX1, SLC2A2 and TCF7L2 SNPs with OS across the populations tested and considering that gender-specific genetic alterations might influence MM survival [48], we suggest that these variants might also exert a modest effect to modulate patient survival.
SCL30A8 gene encodes a zinc transporter involved in the control of insulin processing and secretion [70]. Although no previous studies have reported a link between this locus and MM, there are evidences that suggest that zinc transporters might contribute to cancinogenesis [71,72] through a genderdependent mechanism [73]. The association of the coding SLC30A8 rs13266634 SNP with OS was due to a non-diabetogenic allele suggesting that, rather than modulating glucose homeostasis and insulin secretion [74], the effect attributed to the SLC30A8 locus on MM survival might be driven by a direct effect of Zinc in biological processes such as DNA and RNA stabilization [75], binding of protooncogenes to DNA [75][76][77] and the activation of IGF1 [26,27,61] or telomerase [78]. The SLC30A8 rs13266634C allele has been consistently associated with decreased rates of Zinc transport activitiy and reduced intragranular Zinc levels [79]. However, eQTL data on plasma cells from MM patients did not reveal correlation between this variant and SLC30A8 mRNA levels suggesting that, rather than regulating gene expression, the T allele affect transporter activity in an allele-dose-dependent manner causing increased Zinc concentration and thereby promoting unlimited proliferation of MM cells, disease progression and poor survival. In addition, regulatory data suggest that the SLC30A8 locus might play a role in survival through the modulation of specific transcription factors implicated in tumour promotion and dissemination. www.impactjournals.com/oncotarget As for the HNF1B and SLC30A8 SNPs, the association of the BCL11A rs10490072 and MADD rs7944584 SNPs with OS was determined by non-diabetogenic alleles. BCL11A functions as a myeloid and B-cell protooncogen and has been associated with the development of B-cell malignancies [80,81] whereas MADD encodes for a MAP-kinase activating cell domain involved in the control of physiological cell death through TNFand caspase-dependent apoptosis [82]. In contrast to these associations, the association of the TCF7L2 rs7903146, SLC2A2 rs11920090 PRC1 rs8042680 , and PROX1 rs340874 SNPs with OS was driven by diabetogenic alleles, which suggested that the effect of these variants on OS might be explained by their regulatory effect on insulin secretion and, consequently, on cell proliferation and tumour cell growth. Whereas SLC2A2 encodes a highly efficient glucose transport that is expressed in pancreatic cells and regulates insulin secretion by modulating entry of glucose into the pancreatic cell [83], TCF7L2, PRC1 and PROX1 are proteins that have been involved in β-cell survival and function [84] and in glucose and nonesterified fatty acids or branched-chain amino acids metabolism in liver [85,86]. Despite these interesting results and the regulatory data observed for all these SNPs, the lack of information regarding T2D status among MM patients did not allow us to ensure that the observed effect of the TCF7L2, SLC2A2, PRC1, and PROX1 SNPs on OS could not be due to a different distribution of diabetic patients when grouping by genotype or gender.
This study had both strengths and limitations. Strengths include the use of relative large discovery and replication populations that allowed us to validate the most interesting associations. Limitations include lack of information regarding the classical genetic prognostic factors (chromosomal abnormalities, etc.), T2D status and a relatively small statistical power to detect modest associations with OS, especially when gender-stratified analysis were performed. Another limitation was the partial overlapping of genetic information between studies and the use of imputed genotypes that did not allow to perform a reliable validation of the association observed for genetic variants within ADAMTS9, BCL11A, KCNQ1, MADD and PROX1 genes with overall survival.
In conclusion, this study reports the first evidence of an association between the HNF1B rs7501939 SNP and OS for MM and suggests that the HNF1B locus might, likely through a non-insulin-dependent mechanism, play an important role in modulating MM prognosis. Likewise, this study shows a strong association of the SLC30A8 rs13266634 SNP with poor OS in men that might, at least in part, account for gender differences in OS. Additional studies using larger and well-characterized populations are needed to further replicate these findings but also those involving the BCL11A, MADD, PRC1, PROX1, SLC2A2 and TCF7L2 loci on OS.

Patients, clinical data collection and survival endpoint definition
A total of 1420 Caucasian MM patients were ascertained through the IMMEnSE consortium. Full details of this consortium have been published elsewhere [87]. In brief, inclusion criteria were newly diagnosed MM with Salmon & Durie stage I, II and III, age 18-90 years inclusive and Caucasian origin. DNA was purified from blood specimens using the QIAamp DNA Blood Mini Kit (Qiagen) and clinicophathological characteristics including age, gender, country of origin and disease stage (Durie-Salmon) were retrospectively gathered from medical records in each participant institution ( Table 1). Diagnosis of patients with symptomatic MM was carried out by hematologists according to the International Myeloma Working Group (IMWG) criteria [88,89]. All patients within the IMMENSE consortium for whom survival information was available were included in the study (936 MM cases, 454 women and 482 men) ( Table 1). All participants gave their written informed consent to participate in the study.

SNP selection and genotyping
Fifty-eight variants were selected based on the GWAS for T2D [84, and were genotyped in the IMMEnSE consortium population (Table 4). We considered only SNPs that were replicated in large and independent populations or which came up in several GWAS or their meta-analyses. Additional criteria were potential functionality and linkage disequilibrium (LD) between the reported SNPs. The genotyping of the selected polymorphisms was carried out at GENYO (Centre for Genomics and Oncological Research: Pfizer/University of Granada/Andalusian Regional Government, Granada, Spain) using KASPar ® assays (LGC Genomics, Hoddesdon, UK) according to manufacturer's instructions. For internal quality control, 5% of samples were randomly selected and included as duplicates. Concordance between the original and the duplicate samples for the 58 SNPs was ≥ 99.0% . Call rates for all SNPs were ≥ 90.0% with the exception of the WFS1 rs734312 SNP that was excluded from further analyses.

Replication
For replication purposes, seven hundred MM patients (296 women and 404 men) were provided by the University Clinic of Heidelberg (Germany). This cohort consists of 98 GMMG-HD3 trial patients, 291 GMMG-HD4 trial patients and 311 patients transplanted in Heidelberg but not enrolled in clinical trials (Table 1). Ethical approval for these patients and written informed 1 T/T genotype was associated with poor OS in women with an opposite but not significant effect in men. 2 C allele was associated with better OS in women with no effect in men. 3 G allele was associated with better OS in women with no effect in men. 4 C/C genotype was associated with better OS. No gender-specific effect was observed. 5 T/T genotype was associated with poor OS. No gender-specific effect was observed. 6 T allele was associated with poor OS. No gender-specific effect was observed. 7 T/T genotype was associated with poor OS in women with an opposite but not significant effect in men. 8 C allele was associated with poor OS in men with an opposite but not significant effect in women. 9 T allele was associated with poor OS. No gender-specific effect was observed. 10 T/T genotype was associated with better OS in men with no significant effect in women.

11
G allele was associated with better OS in men with an opposite but not significant effect in women. 12 G allele was associated with poor OS in women with an opposite but not significant effect in men. 13 The presence of each additional copy of the T allele was associated with poor OS in men with no effect in women (additive effect). 14 T allele was associated with better OS in men with no effect in women. 15 T allele was associated with better OS in men with an opposite but not significant effect in women.

eQTL analysis
We also assessed whether selected SNPs correlated with mRNA expression levels in a public eQTL browser for peripheral blood mononuclear cells (http:// genenetwork.nl/bloodeqtlbrowser/) [54]. Expression quantitative trait loci (eQTL) data on malignant plasma cells of 658 patients from the University Clinic of Heidelberg (Germany) were also available for this study. Detailed information on sample collection and clinicopathological characteristics of MM patients as well as technical details of gene expression analysis have been published elsewhere [127].

Statistical analysis
We used chi-square tests to assess Hardy-Weinberg Equilibrium (HWE) for each SNP among IMMEnSE patients. The primary outcome was OS and the endpoint was defined as death from any cause. Survival time was calculated as the time from MM diagnosis (discovery population) or the first stem cell transplantation (replication population) until the occurrence of the study endpoint, censoring at the date of death or the last observed follow-up time. Association with OS defined as hazard ratio (HR) was calculated for each SNP using Cox regression multivariate analysis adjusted for age, gender, country of origin and Durie-Salmon stage (IMMENSE cohort) or for age, gender and clinical trial (Heidelberg cohort). Association estimates were calculated according to dominant, recessive and log-additive models of inheritance with the major allele as reference for regression analyses (Table 4). We also performed gene-gender interaction analyses to determine whether the association between SNPs and MM OS was of similar magnitude in men and women. Survival function was displayed using the Kaplan-Meier method [128] and survival differences across genotypic groups were analysed using the log-rank test.
In order to account for multiple comparisons, we used the Meff/MeffLi method [129], which calculates the effective number of independent genetic markers analysed (N = 54) on the basis of the spectral decomposition (SpD) of matrices of pairwise LD between SNPs (http://neurogenetics. qimrberghofer.edu.au/SNPSpDlite). In addition, we also considered the number of genetic inheritance models tested (dominant, recessive and log-additive). This resulted in a study-wide significance threshold of 0.00031 ([0.05/54]/3) to keep type I error rate at 5% .
Finally, in order to confirm significant associations, a meta-analysis combining genetic data obtained in the IMMENSE population with those extracted from the GWAS conducted in the Heidelberg cohort was also performed following dominant, recessive and additive models of inheritance. The I 2 statistic was used to assess heterogeneity between both studies and the pooled HR was computed using the random-effect model (assuming that between-study variation might depend on chance or random variation and an individual study effect). Randomeffects models are more conservative than fixed-effects models and give rise to wider confidence intervals (CI), which ensures the reliability of the results even though the data come from studies with a relatively different design. All statistics were calculated using SPSS (v. 20) and STATA (v.12) for MAC.