Genetic Analysis Reveals a Significant Contribution of CES1 to Prostate Cancer Progression in Taiwanese Men

The genes that influence prostate cancer progression remain largely unknown. Since the carboxylesterase gene family plays a crucial role in xenobiotic metabolism and lipid/cholesterol homeostasis, we hypothesize that genetic variants in carboxylesterase genes may influence clinical outcomes for prostate cancer patients. A total of 478 (36 genotyped and 442 imputed) single nucleotide polymorphisms (SNPs) in five genes of the carboxylesterase family were assessed in terms of their associations with biochemical recurrence (BCR)-free survival in 643 Taiwanese patients with prostate cancer who underwent radical prostatectomy. The strongest association signal was shown in CES1 (P = 9.64 × 10−4 for genotyped SNP rs8192935 and P = 8.96 × 10−5 for imputed SNP rs8192950). After multiple test correction and adjustment for clinical covariates, CES1 rs8192935 (P = 9.67 × 10−4) and rs8192950 (P = 9.34 × 10−5) remained significant. These SNPs were correlated with CES1 expression levels, which in turn were associated with prostate cancer aggressiveness. Furthermore, our meta-analysis, including eight studies, indicated that a high CES1 expression predicted better outcomes among prostate cancer patients (hazard ratio 0.82, 95% confidence interval 0.70–0.97, P = 0.02). In conclusion, our findings suggest that CES1 rs8192935 and rs8192950 are associated with BCR and that CES1 plays a tumor suppressive role in prostate cancer.


Introduction
Prostate cancer is the most commonly diagnosed cancer and the second most common cause of death among men, with an estimated 191,930 new cases and 33,330 deaths expected worldwide in 2020 [1]. Standardized clinical management approaches, such as radical prostatectomy (RP), radiotherapy, and androgen deprivation therapy, have led to improved outcomes in patients with prostate cancer. However, prognosis remains heterogeneous, suggesting that genetic factors may contribute to treatment response. Genome-wide association studies (GWASs) have successfully identified more than 100 prostate cancer susceptibility loci to date [2,3]. Further functional studies indicate that these risk loci are often located near genes, and regulate genes involved in carcinogenesis, including cell metabolism (JAZF1 and HNF1B) and DNA repair or cell cycle machinery (MYC, TERT, ATM, and CDKN1B) [4]. A pathway-based analysis using known prostate cancer susceptibility loci highlights the antigen presentation pathway and gene network of lipid metabolism, molecular transport, and small molecule biochemistry, which may contribute to prostate cancer development [5]. However, scanning for an association between genetic variants and the prognosis of prostate cancer is still difficult due to the need for a large sample size and long-term follow-up [6,7]. Despite a reasonably large cohort of 24,023 prostate cancer patients with 3513 disease-specific deaths, no evidence of association was observed between genetic variants and prostate cancer survival [7]. GWASs always focus on the most significant genetic variants, which may miss the loci that confer true effects but do not rank at the top. The biological hypothesis-driven approach allows for targeted evaluation and improves the power to detect significant associations. Several functional variants have been reported to be associated with prostate cancer survival by using this approach [8,9].
The carboxylesterase (CES) gene family encodes major liver enzymes, which are responsible for the hydrolysis of various endogenous substrates, including esters, thioesters, amides, carbamates, and xenobiotics, including toxins and drugs [10]. Five human CESs have been identified, and these enzymes share 39-46% of amino acid sequence identity [11]. Although CESs are expressed in most metabolic organs, indicating their protective roles against xenobiotics, they still exhibit different tissue distribution and substrate specificity. CES1 is mainly expressed in the liver and prefers to hydrolyze substrates containing a bulky acyl group and a small alcohol group, whereas CES2 is abundantly expressed in the small intestine and colon and prefers to metabolize esters with a small acyl group and a relatively large alcohol group [12,13]. CESs also appear to participate in the metabolism of fatty acids and cholesterol esters and play a role in the blood-brain barrier system [14], suggesting that the enzymes they encode for serve pivotal physiological functions. Interestingly, CES gene expression has been reported to be downregulated in certain cancer types as the diseases progress [15,16]. The expression of CES genes has also been shown to correlate with chemosensitivity in colorectal cancer [17,18]. A genetic analysis indicated a significant association between the single nucleotide polymorphism (SNP) rs11075646 in the 5' UTR of CES2 and the response rate and time to progression in patients with cancer treated with capecitabine [19]. Therefore, we hypothesized that CES gene polymorphisms may also contribute to the differences in prostate cancer outcomes.
To date, no study has investigated whether CESs could mediate prostate cancer progression. In the present study, we analyzed SNP genotyping data and imputed unobserved SNPs in CES genes to comprehensively assess their impact on disease recurrence in prostate cancer patients who received RP.

Patient Recruitment and Data Collection
This study included 643 Taiwanese patients who underwent RP for localized prostate cancer at three medical centers in Taiwan: Kaohsiung Medical University Hospital, Kaohsiung Veterans General Hospital, and National Taiwan University Hospital, as described previously [20]. The clinicopathological data were obtained from the patients' medical records. Biochemical recurrence (BCR) was defined as two consecutive prostate-specific antigen (PSA) elevation events of 0.2 ng/mL or more [21][22][23][24]. The protocol was approved by the institutional review board of Kaohsiung Medical University Hospital (KMUHIRB-2013132), and each participant provided written informed consent, in accordance with the ethical guidelines.

SNP Selection and Genotyping
We selected 36 haplotype tagging SNPs within the five CES genes and their 10 kb flanking regions with a threshold of a minor allele frequency (MAF) of >0.03, based on the 1000 Genomes data for Han Chinese in Beijing, China and Southern Han Chinese [25]. Genomic DNA was extracted from peripheral blood, and genotyping was conducted using Affymetrix Axiom Genotyping Arrays at the National Centre for Genome Medicine, Taiwan, as described previously [26]. The prediction of the untyped SNPs was performed using Minimac4 with 1000 Genomes Project Phase 3 East Asian reference panels [27,28]. SNPs were filtered by a MAF >0.03 and a Hardy-Weinberg equilibrium >0.001, resulting in 36 SNPs being genotyped and an additional 442 SNPs being imputed.

Statistical Analysis
Analyses were performed using Statistical Package for the Social Sciences software version 19.0.0 (IBM, Armonk, NY, USA). A two-sided P < 0.05 indicated statistical significance; q values were calculated for multiple test correction to report the false discovery rate [38].

Results
This analysis included 643 patients who underwent RP for localized prostate cancer. Their clinical characteristics were presented in Table 1. Two hundred and twenty-eight (35.5%) patients experienced BCR with a median follow-up time of 51 months. Univariate Cox regression indicated that PSA, Gleason score, stage, and surgical margin were significantly associated with BCR (P < 0.001).
We performed a single-locus Cox regression analysis to assess the associations of 36 genotyped SNPs in the five CES genes with BCR. Three SNPs were found to be associated with BCR (P < 0.05 ,  Table S1), of which CES1 rs8192935 remained noteworthy after multiple test correction (q = 0.036). We sought to identify the SNPs better correlated with BCR through imputation, referencing the 1000 Genomes Project (East Asian population). Of the additional 442 SNPs that passed imputation quality control, six SNPs showed superior associations with BCR compared to rs8192935. The strongest signal was shown by rs8192950, which is in linkage disequilibrium with rs8192935 (r 2 = 0.753). The risk of BCR was significantly increased with the number of CES1 rs8192935 G and rs8192950 G alleles (P = 9.64 × 10 −4 and 8.96 × 10 −5 , respectively, Table 2 and Figure 1). Additionally, these two SNPs remained independently and significantly associated with BCR after adjustment for age, PSA, Gleason score, cancer stage, and surgical margin (hazard ratio (HR) 1.43, 95% confidence interval (CI) 1.16-1.76, P = 9.67 × 10 −4 for rs8192935, and HR 1.50, 95% CI 1.24-1.90, P = 9.34 × 10 -5 for rs8192950, Table 2).   We performed a single-locus Cox regression analysis to assess the associations of 36 genotyped SNPs in the five CES genes with BCR. Three SNPs were found to be associated with BCR (P < 0.05, Table S1), of which CES1 rs8192935 remained noteworthy after multiple test correction (q = 0.036). We sought to identify the SNPs better correlated with BCR through imputation, referencing the 1000 Genomes Project (East Asian population). Of the additional 442 SNPs that passed imputation quality control, six SNPs showed superior associations with BCR compared to rs8192935. The strongest signal was shown by rs8192950, which is in linkage disequilibrium with rs8192935 (r 2 = 0.753). The risk of BCR was significantly increased with the number of CES1 rs8192935 G and rs8192950 G alleles (P = 9.64 × 10 -4 and 8.96 × 10 -5 , respectively, Table 2 and Figure 1). Additionally, these two SNPs remained independently and significantly associated with BCR after adjustment for age, PSA, Gleason score, cancer stage, and surgical margin (hazard ratio (HR) 1.43, 95% confidence interval (CI) 1.16-1.76, P = 9.67 × 10 -4 for rs8192935, and HR 1.50, 95% CI 1.24-1.90, P = 9.34 × 10 -5 for rs8192950, Table 2).  To identify the possible effects of these SNPs, functional annotations were extracted from the HaploReg v4.1. Both rs8192935 and rs8192950 have effects on enhancer histone marks and motif alterations and are eQTL SNPs for CES1 (Table S2). The rs8192935 G and rs8192950 G alleles were associated with lower expression levels of CES1 in 322 testis tissues from the GTEx Project (Figure 2A). However, the SNPs are not correlated with expression in prostate tissues, possibly due to the small sample size ( Figure S1). These results indicated that lower CES1 expression would correlate with a poor prognosis in prostate cancer. According to two prostate cancer studies from Taylor and TCGA, lower expression levels of CES1 were associated with prostate cancer, a higher Gleason score, a more advanced stage, and worse survival in patients ( Figure 2B,C). A meta-analysis of eight cohorts of 2064 prostate cancer patients was performed to further evaluate the prognostic significance of CES1. The results showed that higher CES1 expression was significantly related to better prostate cancer prognosis under a random effects model (HR 0.82, 95% CI 0.70-0.97, P = 0.02, Figure 2D).

Discussion
In the present study, we genotyped haplotype-tagging SNPs in CES family genes, following imputation, to fine-map additional SNPs that may be relevant and comprehensively analyze their association with prostate cancer progression. We found that CES1 rs8192935 and rs8192950 might be a prognostic factor for BCR-free survival in patients with prostate cancer. Furthermore, functional studies revealed that these SNPs are in eQTL affecting the expression of CES1 and are subsequently correlated with tumor aggressiveness and prostate cancer prognosis.

Discussion
In the present study, we genotyped haplotype-tagging SNPs in CES family genes, following imputation, to fine-map additional SNPs that may be relevant and comprehensively analyze their association with prostate cancer progression. We found that CES1 rs8192935 and rs8192950 might be a prognostic factor for BCR-free survival in patients with prostate cancer. Furthermore, functional studies revealed that these SNPs are in eQTL affecting the expression of CES1 and are subsequently correlated with tumor aggressiveness and prostate cancer prognosis.
CES enzymes are primarily localized within the endoplasmic reticulum in many tissues and play key roles in both endobiotic metabolism and the activation/detoxification of xenobiotics [10]. CES1, also known as serine esterase 1 or monocyte esterase and cholesterol ester hydrolase, is abundantly produced in the epithelia of metabolic organs including the liver, lungs, and bladder, indicating its protective role against xenobiotics [39]. Interestingly, the inhibition of CES1 in monocytes was shown to diminish their ability to lyse tumor cells [40]. The increased frequency of deficient CES1 enzyme activity has also been reported in non-Hodgkin lymphoma and B-cell chronic lymphocytic leukemia [41,42]. These findings suggest a possible tumor-cell-killing or surveillance function of CES1. Furthermore, CES1 is also a transcriptional target gene of pregnane X receptor (PXR) [43], a key sensor of the body's defense mechanism against xenobiotics. The activation of PXR was found to markedly lower the concentration of circulating androgens, suppress prostate regeneration, and inhibit the growth of human prostate cancer cells [44]. The role of PXR in the homeostasis of androgens may provide clues to the mechanism underlying the observed association between CES1 and prostate cancer progression. In the present study, we found that carriers of the CES1 rs8192935 and rs8192950 G variants had a worse BCR-free survival. Since CES1 is the major enzyme responsible for the hydrolysis of many clinical drugs, CES1 rs8192935 and rs8192950 have been recognized as important pharmacogenetic regulators of treatment outcomes [45,46]. According to the annotation of HaploReg, these two variants may be functional, as they are located at the enhancer region and eQTL of CES1, and are likely to disrupt transcription factor binding motifs in various cells. Consistently, we found that rs8192935 and rs8192950 G alleles were associated with a decrease in the mRNA expression levels of CES1, and lower CES1 expression showed a poorer prognosis for prostate cancer patients. However, these SNPs did not affect CES1 expression in prostate tissue, probably because of the limited prostate samples in the GTEx database. Therefore, further experimental characterization is required to elucidate the function of these SNPs/CES1 in prostate cancer.
This study has several inherent limitations. All the participants in our cohort are Taiwanese, and the findings may not be generalizable to other ethnic groups. Although multiple test correction was performed, the current results still have to be interpreted with caution. In addition, this is a retrospective study with a moderate sample size, and large studies with prospective designs are needed to validate our findings. Finally, no direct biological experiments were conducted to investigate the exact mechanism of action of CES1 rs8192935 and rs8192950 on prostate cancer progression, which should be explored in the future.

Conclusions
Our results suggest that rs8192935 and rs8192950 may reduce the expression of CES1, resulting in a poor prognosis, and could be potential biomarkers of clinical outcome in prostate cancer patients. However, validation in a larger population and further functional studies are needed to identify the tumor suppressive role of CES1 underlying prostate cancer progression.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/2072-6694/12/5/1346/s1. Table S1: Genotyped SNPs and the P values of their association with BCR after RP. Table S2: Regulatory annotation of CES1 rs8192935 and rs8192950. Figure S1: The correlation of rs8192935 and rs8192950 genotypes with CES1 mRNA expression levels in prostate tissues from the GTEx database.