TIMP3 Gene Polymorphisms of -1296 T > C and -915 A > G Increase the Susceptibility to Arsenic-Induced Skin Cancer: A Cohort Study and In Silico Analysis of Mutation Impacts

Long-term exposure to arsenic may induce several human cancers, including non-melanoma skin cancer. The tissue inhibitor of metalloproteinase (TIMP)-3, encoded by the TIMP3 gene, may inhibit tumor growth, invasion, and metastasis of several cancer types. In this study, we aimed to investigate effects of the TIMP3 -1296 T > C (rs9619311) and -915 A > G (rs2234921) single-nucleotide polymorphisms (SNPs) on skin cancer risk in an arsenic-exposed population, and to evaluate the influence of allele-specific changes by an in silico analysis. In total, 1078 study participants were followed up for a median of 15 years for newly diagnosed skin cancer. New cases were identified through linkage to the National Cancer Registry of Taiwan. A Cox regression analysis was used to evaluate the effects of TIMP3 variants. Transcription factor (TF) profiling of binding sites of allele-specific changes in SNPs was conducted using the JASPAR scan tool. We observed borderline associations between TIMP3 genotypes and skin cancer risk. However, when combined with high arsenic exposure levels, the rs9619311 C allele, rs2234921 G allele, or C-G haplotype groups exhibited a greater risk of developing skin cancer compared to the respective common homozygous genotype group. The in silico analysis revealed several TF motifs located at or flanking the two SNP sites. We validated that the C allele of rs9619311 attenuated the binding affinity of BACH2, MEIS2, NFE2L2, and PBX2 to the TIMP3 promoter, and that the G allele of rs2234921 reduced the affinity of E2F8 and RUNX1 to bind to the promoter. Our findings suggest significant modifications of the effect of the association between arsenic exposure and skin cancer risk by the TIMP3 rs9619311 and rs2234921 variants. The predicted TFs and their differential binding affinities to the TIMP3 promoter provide insights into how TIMP3 interacts with arsenic through TFs in skin cancer formation.


Introduction
Arsenic is a ubiquitous element in the Earth's crust. The general population is exposed to arsenic mainly through contaminated groundwater, such as that which is found in Taiwan, West Bengal, Mongolia, and some regions of the United States [1]. Long-term population. Furthermore, to evaluate the impacts of the SNP variants, TF profiling for the binding sites of allele-specific changes was conducted by an in silico analysis.

Baseline Characteristics by TIMP3 Promoter Genotypes
The three studied SNPs of -1296 T > C (rs9619311), -915 A > G (rs2234921), and -899 T > C (rs2234920) present similar allelic frequencies to those of East Asian populations (https://www.ncbi.nlm.nih.gov/snp). The rs2234920 gene locus was not included in the subsequent analysis because of a low minor allelic frequency (Table S1). The other two SNPs had a high degree of linkage disequilibrium (LD), which formed two main components of a haplotype block, T-A and C-G ( Figure S1).
The distributions of the TIMP3 genotype by SNPs are presented in Table 1. In all three genotypes, each SNP in these study participants complied with the Hardy-Weinberg equilibrium (HWE, p > 0.05). As shown in the table, the distribution of baseline characteristics, including age, gender, education level, cigarette smoking, alcohol consumption, and arsenic exposure, were similar among the three genotypes for each respective SNP (p > 0.05 for all).

Factors Associated with the Skin Cancer Incidence
In total, 50 patients were newly diagnosed with skin cancer during a median 15-year follow-up. As shown in Table 2, an older age, male gender, and higher levels of arsenic exposure were significantly associated with an increased risk of skin cancer in the study participants. On the other hand, a higher education level and habitual smoking were associated with a decreased risk of skin cancer. Table 3 shows the relationship between the TIMP3 promoter genotype and the risk of skin cancer. Although not statistically significant, there was a trend toward a positive association between the skin cancer risk with carriers of the SNPs of the rs9619311 C allele and rs2234921 G allele. In the haplotype analysis, no single diplotype was significantly associated with the risk of skin cancer (Table S2).

Synergetic Effect of the TIMP3 Genotype and Arsenic Exposure
Although the TIMP3 genotype was a weak independent predictor of skin cancer risk, we examined whether the genotype and arsenic exposure had an interactive effect on the skin cancer risk. As shown in Table 4, when combined with high levels of arsenic exposure, the rs9619311 C allele (HR = 3.29, 95% CI = 1.07~10.12, dominant model), rs2234921 G allele (HR = 3.81, 95% CI = 1.30~11.23, dominant model; HR = 82.52, 95% CI = 8.60~791.61, recessive model), and C-G haplotype (HR = 3.31, 95% CI = 1.08~10.18, dominant model) groups had a greater risk of developing skin cancer compared to the respective common homozygous genotype groups.

TIMP3 Genotype and Incidence of Skin Lesions
We further examined the relationship between TIMP3 promoter genotypes and the incidence of skin lesions in the study participants. There was an increased risk of hyperpigmentation (OR = 6.72, 95% CI = 1.59~28.42) in the rs2234921 G/G genotype group, as shown in Table 5. The G allele of rs2234921 of the TIMP3 gene was only marginally associated with a hyperkeratosis risk (p = 0.06). We also found a significant association between arsenic exposure and the risk of skin lesions in the rs2234921 G allele group (HR = 2.96, 95% CI = 1.25~6.96, dominant model; HR = 19.82, 95% CI = 4.34~90.57, recessive model), who had a greater risk of developing skin lesions compared to the common homozygous genotype group (Table 6). These results also support a synergetic effect of the TIMP3 genotype and arsenic exposure on the skin lesion risk. * All variables are categorized and presented as counts (percentage). † Differences between total counts and the total number (n) by genotype are due to missing data.

In Silico Analysis of TF-Binding Sites (TFBSs) and the Impacts of Mutations
To discover potential TF motifs located at or flanking the rs9619311 and rs2234921 SNPs, we launched a detailed survey using a eukaryotic promoter database [23] and the JASPAR database [24]. Among 568 TF motifs, rs9619311 and rs2234921 were found to respectively be located within the predicted sequence of 15 and 3 transcription motifs ( Figure S2). To further analyze whether the T to C conversion of rs9619311 affected the binding affinities of each motif, we used the JASPAR scan tool to obtain TFBS scores and compared the scores between the wild-type (T) and mutant (C) rs9619311. Accordingly, mutant rs9619311 was found to have attenuated binding to the BACH2, MEIS2, NFE2L2, and PBX2 TFs as shown in Figure 1. Mutant rs2234921 (G allele) was also found to have a reduced binding affinity of E2F8 and RUNX1 to their motifs, as shown. Interestingly, we did not find differences in the binding affinities of NF-κB by rs9619311, although a putative NF-κB motif flanking rs9619311 was identified. In addition, we found no AP-1 motif in rs9619311 or rs2234921.    [25], MEIS2 [26], NFE2L2 [27], PBX2 [28], E2F8 [29], RUNX1 [30]. N.A., not available.

Discussion
In this study, we examined contributions of the TIMP3 -1296 T > C (rs9619311) and -915 A > G (rs2234921) SNPs to arsenic-induced skin cancer/lesion risk. In the recessive model, we observed marginal and significant associations of the G allele of the rs2234921 SNP with the risk of arsenic-induced skin cancer and skin lesions, respectively. The rs9619311 SNP alone or in the haplotype block showed no association with the skin cancer or skin lesion risk. Although there were few associations of the studied SNPs with the skin cancer/lesion risk, we found a significant modification effect of variant alleles on the association between arsenic exposure and skin cancer/lesion risk in either the dominant or recessive model.
Previous studies on the associations between the TIMP3 -1296 T > C and the risk of various cancers produced inconsistent results. A case-control study reported that variant genotypes (T/C or C/C) were positively associated with breast cancer risk [17]. In contrast, other studies reported associations of variant genotypes with a decreased risk of cancer, such as hepatocellular carcinoma in females and colorectal cancer [15]. However, the most recent case-control studies or patient survival studies have so far found no associations between variant genotypes and the risk of cancers, including breast [18], bladder [21], gastroesophageal junction [16], and oral cavity cancers [19]. Our study did not find a significant association between the TIMP3 -1296 T > C SNP and skin cancer or skin lesion risks, although a borderline association of the -915 A > G SNP, an SNP in LD with the -1296 T > C site, was observed. To our knowledge, no observational studies have investigated the associations of the -915 A > G SNP with various cancers. Discrepancies in the literature may be attributed to different types of cancers, different stages of the same cancer, or other yet unidentified factors.
In the present study, we found a significant effect modification of the studied SNPs on the association between arsenic exposure and skin cancer/lesion risk. Subjects who carried a variant genotype (T/C or C/C of rs9619311 or A/G or G/G of rs2234921) were at an increased risk of skin cancer/lesions that were associated with arsenic exposure compared to those who carried the respective common homozygous genotypes (T/T of rs9619311 and A/A of rs2234921). To our knowledge, this study is the first report of an interaction between TIMP3 genetic polymorphisms and an environmental contaminant such as arsenic. These results suggest that genetic variants of the TIMP3 promoter region could play a role in arsenic-induced skin carcinogenesis.
Arsenic is a strong oxidant that rapidly induces several stress proteins in cultured cells [31]. During arsenic metabolism in cells, reactive oxygen species are generated [32]. Although the exact molecular mechanisms of arsenic carcinogenesis are not well understood, it is generally accepted that free radicals that are generated by arsenic may induce intracellular signal transduction and activate redox-sensitive transcription factors such as AP-1, NF-κB, and nuclear factor erythroid 2-related factor 2 (NRF2), which in turn change expressions of genes that are involved in cell growth, proliferation, and malignant transformation [33,34]. However, in this study, the in silico analysis revealed that the two studied SNPs did not predict the binding sites for either AP-1 or NF-κB. This finding is consistent with a previous report by Beranek et al. [22]. Interestingly, we found several TF motifs located at or flanking the two SNP sites of TIMP3 that had not previously been reported. We further validated that the C allele of rs9619311 attenuated the binding affinity of BACH2, MEIS2, NFE2L2, and PBX2 to the TIMP3 promoter. We also validated that the G allele of rs2234921 reduced the affinity of E2F8 and RUNX1 to bind to their site on the promoter. BACH2 was identified as a competitor for the Maf protein in the NRF2 antioxidant defense pathway [35] and was recently shown to be highly sensitive to DNA damage and aging [25]. MEIS2 is a homeobox protein usually in complex with PBX and HOX to promote downstream gene transcription that is involved in early development and cell differentiation [26]. Others also showed a correlation between MEIS2 and the progression of various cancers [36][37][38][39]. NFE2L2 encodes NRF2, which is the master regulator of cell defense against stress [27], suggesting that TIMP3 may play a role in downstream gene responses to cellular redox imbalances. PBX2 is a member of the PBX family that was initially identified as an essential Hox cofactor that plays crucial roles in early development [28]. Recent studies suggested that PBX2 dysregulation is closely associated with cancer progression [40,41]. E2F8 is a transcription repressor that antagonizes E2F1 in regulating the cell cycle, apoptosis, and tumor promotion [29]. RUNX1 was initially described as a critical regulator of developmental hematopoiesis [30]. Recently, RUNX1 was also associated with solid tumor development in the skin, lungs, intestines, and breast [42]. Taken together, the reduced affinity of the above TFs to TIMP3 promoter variants may result in decreased TIMP3 gene expression and activity. Although the predicted TFs and their differential binding affinities need to be confirmed by experimental approaches, our results agree with the reported functional effects of TIMP-3 on essential cellular processes, including proliferation, differentiation, and apoptosis in the regulation of cancer genesis.
The strength of this study is the prospective study design with newly diagnosed cancer/lesion cases. Unlike the case-control study design of most previous TIMP3 SNP studies, the possibility of reverse causation could thus be minimized. In addition, this study used a national registry database to determine the incident cancer cases, which contains accurate and complete datasets that are periodically updated. Several limitations of this study should be noted. First, because of the small sample size, we could not evaluate the cancer risk by cancer subtype, such as Bowen's disease vs. invasive skin cancer. Whether our findings differ among subtypes is not known. The small number of incident cases is also of concern; hence this study's findings may have been due to chance. In addition, most of our results showed a modest to large effect size of the association but did not reach statistical significance. One explanation may be due to the lower minor allele frequency (MAF) observed in the two SNPs examined in our study participants. However, the MAF occurred at a similar frequency to that reported in previous studies involving Taiwanese [19,20] and as reported in the dbSNP database that includes East Asians (https://www.ncbi.nlm.nih.gov/snp). The allelic and genotypic frequencies of our study participants should have been well represented.

Study Participants and Questionnaire Data
The recruitment of the study participants and the collection of questionnaire data and blood samples were previously described [5,9,43]. In brief, study participants were recruited from two arseniasis endemic areas in Taiwan: a Blackfoot disease-endemic area in the southwestern region (LMN subcohort) and the Lanyang Basin in the northeastern coastal region (Lanyang subcohort). Residents of these two areas were exposed to arsenictainted groundwater in the 1910s~1970s and 1940s~1990s, respectively. Arsenic exposure is associated with increased risks of non-melanoma skin cancer [44,45], lung cancer [43], and urinary tract cancer [46]. Epidemiologic follow-up studies were launched in 1989 and 1998 to identify susceptibility factors for arsenic-associated diseases. Since then, periodic physical examinations to identify skin cancer and skin lesions have been carried out in the endemic areas [4,5,47].
This study consisted of 1078 study participants from the two study subcohorts who had baseline questionnaire data and successful genotyping results. These study participants were reported to have increased mortality from cardiovascular diseases and an increased incidence of a variety of cancers associated with arsenic [9,48]. A prior study identified demographic and lifestyle risk factors for cancers in the endemic areas, such as age, gender, education level, cigarette smoking, and alcohol consumption which were considered to be potential confounding factors in the present study for the association analysis [9]. Individuals were considered regular users if they smoked cigarettes or consumed alcohol at least 3 days a week for at least 6 months. To represent the overall exposure to arsenic from artesian well water for each study subject, cumulative arsenic exposure was applied as an index of arsenic exposure, calculated as previously described [4]. Informed consent was obtained for participation at the time of enrollment [9,49]. This study was performed in accordance with the Declaration of Helsinki and was approved by the institutional review boards of Taipei Medical University (N201807031, 2018 July 24) and the Genomics Research Center, Academia Sinica (AS-IRB01-09008, 6 April 2009 and AS-IRB01-11070, 3 March 2012).

Identification of Skin Cancer Cases
We identified newly diagnosed skin cancer patients through a data linkage with profiles of the National Cancer Registry of Taiwan and Death Certification System of Taiwan as previously described [9]. The nationwide Cancer Registry was implemented in 1979. Both invasive cancers and carcinoma in situ of the skin (Bowen's disease) were included in the present study. In Taiwan, it is mandatory to report Bowen's disease and cervical carcinoma in situ to the Cancer Registry. The percentage of pathological confirmation of skin cancer cases was 99.7% [50]. The Death Certification System has been computerized since 1972, and the system maintains updated information on the vital status of Taiwanese citizens. The total number of person-years (P-Y) in the follow-up period for each participant was calculated from the interview date to the date of cancer diagnosis, death, or the end of follow-up (31 December 2014), whichever occurred earliest.

Identification of Skin Lesion Cases
Skin lesions, including hyperpigmentation and hyperkeratosis, were clinically diagnosed during community health examinations held from 1989 to 2009 for the LMN subcohort. The total number of person-years for each participant was calculated from the date of the skin examination with no manifestation of skin lesions to the date of lesions being newly diagnosed, death, or the end of follow-up (30 April 2009), whichever occurred earliest. The vital status was determined by linking it to the Death Certification System as described above.

TIMP3 Promoter Genotyping
DNA were extracted from the buffy coat according to the manufacturer's instructions (Gentra System, Minneapolis, MN, USA). The primer sequences that were used for fragment amplification were 5 -GAAAGGGGTGACGAGTTCCTG-3 (forward) and 5 -CCTTGACTGTGCTTGGTGGA-3 (reverse). The polymerase chain reaction (PCR) analysis was carried out in a final volume of 25 µL containing 12.5 µL of 2× Tag PCR mix (TIANGEN Biotech, Beijing, China), 1 µL of 10 µM of each primer, 0.5 µL of 10 mM dNTP, and 1 µL of a 5-ng DNA sample. The thermal cycling conditions were as follows: pre-denaturation at 94 • C for 5 min, 33 cycles of denaturation at 94 • C for 30 s, annealing at 65 • C for 30 s, and extension at 72 • C for 30 s. Finally, the process was extended at 72 • C for 5 min. All the PCR products (903 bp) were sent to the Institute of Biomedical Sciences, Academia Sinica (Taipei, Taiwan), for a sequencing analysis ( Figure S3). A total of five percent of the analyzed samples were randomly selected for a reproducibility analysis, and the repeated results were 100% consistent.

In Silico Analysis of TF-Binding Profiles
The EPD server (https://epd.epfl.ch//index.php, accessed on 5 October 2022) is an in silico tool for analyzing promoter-binding regions [23]. The Search Motif Tool was applied to analyze TF motifs located in TIMP3 based on bioinformatics from JASPAR (https://jaspar.genereg.net, accessed on 5 October 2022). JASPAR is an online server for profiling TF binding, position frequency matrices (PFMs), and TF flexible models (TFFMs). The score of each binding site was calculated as log(P(Site|Matrix)/P(Site|Background)) [24].

Statistical Analysis
For the genetic analysis, we first examined the Hardy-Weinberg equilibrium (HWE) in the study participants using the χ 2 goodness-of-fit test. We then used Haploview v.4.2 (Broad Institute, Cambridge, MA, USA), to measure linkage disequilibrium (LD) between the studied SNPs and to build a haplotypic block. To compare the frequency distributions of baseline characteristics among the genotype groups, we used the χ 2 test or Fisher's exact test. To analyze the associations of the studied SNPs with skin cancer/lesion risks, we used Cox proportional hazard models while adjusting for other covariates. All variables were categorized except for age. Arsenic exposure was categorized into three groups ≤ 7.7, 7.7~17.5, and > 17.5 ppm-years with reference to previous reports [4,51]. Trend tests were estimated by treating the three exposure groups as a continuous variable, and integer scores were coded in the regression analysis. To examine the association of TIMP3 genotypes with skin cancer/lesion risk, we calculated hazard ratios (HRs) and 95% confidence intervals (CIs) derived from Cox regression analyses based on three genetic models: additive, dominant, and recessive.
To assess whether there was an interaction between the TIMP3 promoter genotype and arsenic exposure on skin cancer/lesion risk, we conducted a combination analysis of the two variables on the cancer/lesion risk. We examined whether the group with the risk genotype and high arsenic exposure (>17.5 ppm-years) exhibited a higher risk than the group with the risk genotype alone or high arsenic exposure alone. Additive interactions and multiplicative interactions were respectively evaluated using a trend test and a likelihood ratio test. All statistical analyses were performed using SAS 9.4 (SAS Institute, Cary, NC, USA), and two-tailed p < 0.05 was considered significant.

Conclusions
Based on this cohort study, our data suggest that TIMP3 -1296 T > C and -915 A > G promoter SNPs are not likely significant predictors of skin cancer/lesion risk among Taiwanese participants. However, we found a significant effect modification of the association between arsenic exposure and skin cancer/lesion risk by variant alleles of the studied SNPs. These findings have implications for risk stratification in arsenic risk assessments. Our results of the in silico analysis provide insights into the roles of several transcription factors that bind to the two TIMP3 promoter SNP sites. These findings may facilitate further investigation of how TIMP3 interacts with arsenic through transcription factors and its role in cancer formation. Informed Consent Statement: Informed consent was obtained from all subjects that were involved in the study. Written informed consent was obtained from the subjects to publish this paper.
Data Availability Statement: Datasets of the National Cancer Registry and Death Registry System analyzed and/or generated during this study are held by the Ministry of Health and Welfare (MOHW) of Taiwan. The data are available from the corresponding author with the permission of the Taiwan MOHW.