Rs10204525 Polymorphism of the Programmed Death (PD-1) Gene Is Associated with Increased Risk in a Saudi Arabian Population with Colorectal Cancer

Checkpoint programmed death-1 (PD-1) has been identified as an immunosuppressive molecule implicated in the immune evasion of transformed cells. It is highly expressed in tumor cells in order to evade host immunosurveillance. In this study, we aimed to assess the association between single nucleotide polymorphisms (SNP) of PD-1 and the risk of colorectal cancer (CRC) in the Saudi population. For this case-control study, the TaqMan assay method was used for genotyping three SNPs in the PD-1 gene in 100 CRC patients and 100 healthy controls. Associations were estimated using odds ratios (ORs) and 95% confidence intervals (95% CIs) for multiple inheritance models (codominant, dominant, recessive, over-dominant, and log-additive). Moreover, PD-1 gene expression levels were evaluated using quantitative real-time PCR in colon cancer tissue and adjacent colon tissues. We found that the PD-1 rs10204525 A allele was associated with an increased risk of developing CRC (OR = 2.35; p = 0.00657). In addition, the PD-1 rs10204525 AA homozygote genotype was associated with a high risk of developing CRC in the codominant (OR = 21.65; p = 0.0014), recessive (OR = 10.97; p = 0.0015), and additive (OR = 1.98; p = 0.012) models. A weak protective effect was found for the rs2227981 GG genotype (OR = 2.52; p = 0.034), and no significant association was found between the rs2227982 and CRC. Haplotype analysis showed that the rs10204525, rs2227981, rs2227982 A-A-G haplotype was associated with a significantly increased risk of CRC (OR = 6.79; p =0.031).


Introduction
Colorectal cancer (CRC) is a malignant neoplasm that arises from the lining of the large intestine. It appears following the accumulation of multiple genetic and epigenetic alterations in cells, leading to many abnormalities in several central signaling pathways [1,2]. Both exogenous and endogenous factors contribute to the increased risk of CRC. Environmental factors including a low fiber diet, low physical activity, obesity, smoking, and alcohol consumption are all considered to be risk factors for CRC [3][4][5]. Moreover, inherited genetic factors have been shown to contribute significantly to the susceptibility of CRC, as for many other cancer diseases [6][7][8][9]. In addition, the risk of CRC is higher in individuals with inflammatory bowel diseases such as Crohn's disease, which involve alterations in immune responses [10].
Genome-wide association studies (GWASs) have improved the search for genetic variations such as single nucleotide polymorphisms (SNPs) which are related to many diseases including cancer. A significant number of SNPs have been reported as a part of disease susceptibility to CRC. In the European population, GWASs have reported almost

SNP Selection and Genotyping
For this genetic association study, three PD-1 SNPs were selected; rs10204525 (PD-1.6), rs2227981 (PD-1.5), and rs2227982 (PD-1.9). SNPs were acquired from Thermo Fisher Scientific (Thermo Fisher Scientific, Waltham, MA, USA). The characteristics of each SNP are presented in Table 2. Genotyping was performed using a TaqMan allelic discrimination assay as previously described [38]. Genotyping were performed according to the manufacturer's protocol. Briefly, for each PCR reaction, a 20-ng DNA sample was used with 10 µL of 2X Universal Master Mix and 1X assay mix in a total reaction volume of 20 µL (Applied Biosystems, Foster City, CA, USA). PCR conditions were as follows: pre-read stage at 60 • C for 30 s, hold stage at 95 • C for 10 min, PCR stage at 95 • C for 15 s and 60 • C for 1 min for 40 cycles, and post-read stage at 60 • C for 30 s. All genotypes were determined using end-point reading on a ViiATM 7 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). For quality control, 5% of samples were randomly selected and subjected to repeat analysis as a measure of verification for the genotyping procedures. The results were reproducible without any discrepancies. For each DNA sample, allelic genotyping was used for the detection of three SNPs in the PD-1 gene. Relative quantification (RQ) was calculated using the comparative CT method (2(−∆∆Ct)) [38]. Normal tissue samples were used as a calibrator, and GAPDH as a reference gene for normalization. The data were expressed as median values.

Statistical Analysis
Association analysis including the codominant, dominant, recessive, over-dominant, and log-additive models were performed using web-based software (SNPStats) [44]. The inheritance model with the lowest Akaike information criteria (AIC) value was considered as the best fit. Odds ratios (OR) with 95% confidence intervals (CI) were calculated based on logistic regression. Haplotype analysis and linkage disequilibrium (LD) were also performed using SNPStats and Haploview software (http://www.broadinstitute.org/ haploview/haploview/index.php, accessed on 29 August 2022).
Haplotype frequency was estimated through a statistical method to assess the differences in haplotype frequency distribution between cases and controls. For each SNP, deviation from a Hardy-Weinberg equilibrium and χ 2 values were calculated using the web-based tool at https://ihg.helmholtz-muenchen.de/cgi-bin/hw/hwa1.pl, accessed on 2 September 2022.
Haldane-Anscombe correction was applied to estimate the OR for zero values. A p-value of <0.05 was considered as significant.

Demographic Characteristics of Study Population
Patient baseline characteristics are shown in Table 1. The study included 200 participants, 100 patients with sporadic colorectal cancer and 100 healthy matched controls. The patient group comprised 64 males and 36 females with an average age group of 56.33 ± 14.56. The control group comprised 65 males and 35 females with an average age group of 56.31 ± 14.56 (age-and gender-matched controls). Patients were classified according to the three TNM stages: I, II, and III. In total, 57% belonged to stage II and 32% to stage III, and only 11% to stage I. For all subjects, genotyping of the PD-1 gene was performed using a TaqMan assay for three selected SNPs; PD-1.6, PD-1.5, and PD-1.9 (Table 2).

PD-1 SNP Association with CRC
The genetic and allelic association between the three SNPs and CRC was tested using five genetic models (allelic, codominant, dominant, recessive, over-dominant, and additive). The results are presented in Table 3. The distribution of genotypes for PD-1.6 and PD-1.5 in the control group followed a Hardy-Weinberg equilibrium (p > 0.05, χ 2 < 3.84), while the PD-1.9 (p = 0.016, χ 2 = 8.22) was deviated from a Hardy-Weinberg equilibrium.
The genotyping of PD-1.6 showed that the AA genotype was present in 9% of the patients; however, it was lacking among controls. Statistically, the AA genotype was significantly associated with an increased risk of CRC in the codominant (OR = 21.65; 95% CI (1.23-378.73); p = 0.0014), recessive (OR = 10.97; 95% CI (1.37-87.43); p = 0.0015), and additive (OR = 1.98; 95% CI (1.14-3.45); p = 0.012) models, after applying Haldane-Anscombe correction. The A allele was significantly high in frequency in CRC patients compared to the healthy control group (0.19 vs. 0.09), suggesting an increased susceptibility to CRC for individuals sharing this allele (OR = 2.35; 95% CI (1.24-4.03); p = 0.00657). For the PD-1.5 polymorphism, our analysis showed that the AA genotype was the most frequent in both patients and controls, followed by the AG and GG genotypes, respectively. In the codominant model, the GG genotype was found with a low frequency in CRC (8%) compared to healthy controls (18%), although this difference did not reach significance (OR = 2.25; 95% CI (0.90-5.63); p = 0.078). However, in the recessive model, a positive correlation with disease was observed for AG + AA vs. G/G (OR = 2.52; 95% CI (1.04-6.11); p = 0.034) with the lowest AIC (276). Furthermore, the A allele was the most common allele in comparison to the G allele, which does not coincide with the MAF database (A = 0.33/72 (Qatari), A = 0.35/1759 (1000Genomes) ( Table 1). The dominant genotype AA against the AG + GG genotype had a frequency of 52 to 48 with no differences between patients and controls.
Our analysis for PD-1.9 did not show any significant association with CRC in any of the studied models.

Age and Gender Stratified Analysis
In order to investigate the possible influence of the SNPs according to age and gender, we have stratified the 100 cases into two subgroups for each parameter. For age-stratified analysis, we classified the patients to those aged ≥ 56 (n = 59) and those aged < 56 years (n = 41). For gender stratification, the group of females comprised 36 individuals versus 64 males. The allelic and genetic distribution of the three SNPs among each group was evaluated using the five models of inheritance. Our analysis revealed no significant associations between the three SNPs in relation to age (Table 4) or gender (Table 5).

Haplotype Analysis
Linkage disequilibrium (LD) analysis in the control samples revealed a weak LD between PD-1.5 and PD-1.9, and no LD for PD-1.6 ( Figure 1). The haplotypes were generated using the three SNPs among cases and controls (Table 6). Six common haplotypes of PD-1.6, PD-1.5, and PD-1.9 (frequency > 1.8%) showed an accumulated frequency of more than 95% of haplotypes. The distribution of haplotypes exhibited differences between cases and controls. The most frequent haplotype was the G-A-G, with 28% in CRC cases and 41% in healthy controls, which was used as a reference. The haplotype A-A-G was associated with an increased risk of CRC (OR = 6.79; 95% CI (1.21-38.25); p = 0.031). Globally, an association between haplotypes and diseases was supported by the global p-value (p-value = 0.017) ( Table 6).
of PD-1.6, PD-1.5, and PD-1.9 (frequency > 1.8%) showed an accumulated frequency of more than 95% of haplotypes. The distribution of haplotypes exhibited differences between cases and controls. The most frequent haplotype was the G-A-G, with 28% in CRC cases and 41% in healthy controls, which was used as a reference. The haplotype A-A-G was associated with an increased risk of CRC (OR = 6.79; 95% CI (1.21-38.25); p = 0.031). Globally, an association between haplotypes and diseases was supported by the global p-value (p-value = 0.017) ( Table 6).

Gene Expression Analysis of PD-1 mRNA
The quantification of PD-1 gene expression from 35 colon cancer fresh tissues and 35 normal adjacent matching tissues was performed using real-time reverse transcriptase PCR (qRT-PCR). The expression of PD-1 mRNA was significantly higher (4-8-fold difference; p < 0.001) in colon cancer tissues as compared to healthy adjacent colon samples (Figure 2).

In Silico Functional Analysis of the PD-1.6 Polymorphism
In order to evaluate the functional effect of the rs10204525 polymorphism located in the 3 UTR of the PD-1 gene, we have performed in silico analysis. Using the MicroSNiPer prediction website, we have identified the three most plausible miRNAs which possess a 8-10-nt seed length and recognize the miRNA response element (MRE), enclosing the PD-1.6 site ( Table 7). The MFE and ∆MFE are the lowest for hsa-miR-541-3p and hsa-miR-4717-3p, and thus these could be considered as potential miRNAs that regulate gene expression through the MRE that includes the PD-1.6 polymorphism (Figure 3). A genomewide annotation of variant analysis (GWAVA) predicted a deleterious effect of the PD-1.6 mutation with 78% accuracy.

Gene Expression Analysis of PD-1 mRNA
The quantification of PD-1 gene expression from 35 colon cancer fresh tissues and 35 normal adjacent matching tissues was performed using real-time reverse transcriptase PCR (qRT-PCR). The expression of PD-1 mRNA was significantly higher (4-8-fold difference; p < 0.001) in colon cancer tissues as compared to healthy adjacent colon samples ( Figure 2). Figure 2. PDCD1 gene expression assessed by qRT-PCR in colon cancer tissues compared to normal matching tissues: (mean ± SD, data normalized to GAPDH, *** p < 0.00 t-test).

In Silico Functional Analysis of the PD-1.6 Polymorphism
In order to evaluate the functional effect of the rs10204525 polymorphism located in the 3′ UTR of the PD-1 gene, we have performed in silico analysis. Using the MicroSNiPer prediction website, we have identified the three most plausible miRNAs which possess a 8-10nt seed length and recognize the miRNA response element (MRE), enclosing the PD-1.6 site ( Table 7). The MFE and ΔMFE are the lowest for hsa-miR-541-3p and hsa-miR-4717-3p, and thus these could be considered as potential miRNAs that regulate gene expression through the MRE that includes the PD-1.6 polymorphism (Figure 3). A genome-wide annotation of variant analysis (GWAVA) predicted a deleterious effect of the PD-1.6 mutation with 78% accuracy. Table 7. Putative miRNA targeting PD-1.6 SNPs in the 3′ UTR of the PD-1 gene, and the difference in the free energy of hybridization (ΔFME) for wild-type and variant alleles.

Discussion
Functional polymorphisms in immune checkpoint genes are reported to have rious impact on the outcomes of many types of cancers [17,45]. In fact, there is a

Discussion
Functional polymorphisms in immune checkpoint genes are reported to have a serious impact on the outcomes of many types of cancers [17,45]. In fact, there is an increasing focus on the role of SNPs in the modulation of individual susceptibility or protection against cancer. Large-scale GWASs have identified many loci associated with the risk of CRC [6][7][8][9]46]. In this context, certain polymorphisms in the PD-1 gene have been studied and found to contribute to the individual risk of cancer [28,47,48]. Studies have also reported changes in the expression of PD-1 and an association with many types of cancers [21,[49][50][51]. In the current study, we investigated the relationship between the PD-1.6 (rs10204525), PD-1.5 (rs2227981), and PD-1.9 (rs2227982) polymorphisms and the risk of developing CRC in Saudi Arabia. In our study, a strong association between the A allele and AA genotype of PD-1.6 and CRC was found in almost all tested inheritance models. For the two other PD-1 polymorphisms (PD-1.5 and PD-1.9), an association with the disease was found only for the GG genotype of PD-1.5 in the recessive model. For the PD-1 gene, few studies have reported an association between the SNPs tested herein and cancer diseases. The PD-1.6 polymorphism is reported to be associated with an increased susceptibility of some clinical features of esophageal cancer in the Chinese Han population [52]. Additionally, this SNP was correlated with the susceptibility to human T-cell leukemia virus type 1 and some clinical features in an Iranian population [53]. However, in a recent study reported by Fathi, et al. [54], no association between PD-1.6 and basal cell carcinoma (BCC) was found. Similarly, in two other separate case-control studies, no associations were detected between PD-1.6 polymorphism and the risk of developing CRC in subjects from Heilongjiang Province in Northeast China [2] and Eastern China [36]. However, associations between PD-1.6 polymorphisms and several infectious and autoimmune diseases have been reported in many studies [20,27,29,30,32,55,56]. For PD-1.5 and PD-1.9, controversial results have been reported for their association with cancer diseases. Two separate studies reported no association between PD-1.9 and CRC in two populations in China, conducted by Ge et al. [2] and Lin et al. [36]. However, an association between PD-1.5 and an increased susceptibility to colon cancer and rectal cancer was reported earlier by Mojtahedi et al. [37] in an Iranian population. In a study reported by Zhao et al. [22], who examined the association of the three SNPs (PD-1.6, PD-1.5, and PD-1.9) with CRC among a Han Chinese population, no significant connection was found with the disease. In other case-control studies regarding the relationship between PD-1.5 and PD-1.9 and other cancer diseases, including lung adenocarcinoma [31], cervical cancer [32], breast cancer [28], gastric cancer [34], and thyroid cancer [35], significant associations have been reported. This discrepancy in results could be explained by ethnic factors alongside other environmental factors that are usually not considered in the association analysis, such as smoking, diet, alcohol drinking, ultraviolet and ionizing radiation exposition, medical procedures and drugs, occupation, reproductive behavior, pollution, infection, and other still unknown factors [3][4][5]. These environmental factors could act by inducing somatic mutations or epigenetic effects that could modify the structure of DNA or affect its expression or stability [57,58]. In this study, we tested the relationships between the PD-1.6, PD-1.5, and PD-1.9 haplotypes and CRC, and we found a strong positive association between the A-A-G haplotype and an increased risk of CRC. In a study reported by Zhao et al. [22], a protective effect of the A-G-A of the same SNPs against CRC was found among the Han Chinese population. This haplotype was found with a very low frequency in our population. Moreover, we have explored the gene expression levels in colon cancer tissues in comparison to surrounding normal colon tissues. Our results showed an up-regulation of the PD-1 gene in colorectal tumor cells. This result is in agreement with other studies on the association of high-level PD-1 expression with the clinical features of multiple tumors including CRC, hepatocellular carcinoma, esophageal cancer, and kidney clear cell carcinoma [21,26,[49][50][51]59]. Thus, immunotherapy based on monoclonal antibodies against PD-1 proteins, blocking the binding with PD-L1, have demonstrated the efficacy against some malignant diseases [60,61]. Moreover, the PD-1.6 polymorphism has been associated with an overexpression of the PD-1 gene [59]. It has been shown that this polymorphism is involved in the activation and transcription of PD-1. In our in silico analysis, we have identified three putative and allele-sensitive miRNA molecules that could interact with the region enclosing the PD-1.6 polymorphism. miR-541 and miR-4717 were found to be the most plausible; they could play important roles in the regulation of PD-1 gene expression and are sensitive to PD-1.6 polymorphism. In an experimental study, Zhang, et al. [62] showed that miR-4717 may interact in an allele-specific manner with the 3 UTR of PD-1.6 and regulate PD-1 expression. Therefore, lymphocytes from patients with chronic HBV harboring the PD-1.6 GG genotype and transfected with miR-4717 mimics exhibited significantly decreased PD-1 expression and increased TNF-α and IFN-γ production, although no effect was observed with the PD-1.6 AA genotype carriers [62]. Further studies are necessary to investigate the role of PD-1.6 polymorphism and its association with PD-1 expression for eventual use in specific tailored immunotherapy against cancer diseases.

Conclusions
In conclusion, our findings support an association between the PD-1 rs10204525 polymorphism and an increased CRC risk in the Saudi population. Additional reports involving larger sample sizes with more detailed clinical information in different ethnicities will be important to confirm our conclusions.