TGF-β/VEGF-A Genetic Variants Interplay in Genetic Susceptibility to Non-Melanocytic Skin Cancer

Differential genetically determined expression of transforming growth factor-β (TGF-β pathway and of vascular endothelial growth factor-A (VEGF-A) might modulate the molecular “milieu” involved in the etio-pathogenesis of non-melanoma skin cancer (NMSC). We have evaluated the frequency of some functionally relevant SNPs of TGF-β and VEGF-A genes in 70 NMSC patients and 161 healthy controls, typed for TGF-β1 rs1800471, TGF-β2 rs900, TGF-βR1 rs334348 and rs334349, TGF-βR2 rs4522809 and VEGF-A rs3025039 SNPs. TGF-βR2 rs1800629G allele and related genotypes were found to be associated with a possible protective role against NMSC, whereas VEGF-A rs3025039T was associated with an increased risk. To evaluate the effect of genotype combinations on NMSC susceptibility, we determined the frequencies of 31 pseudo-haplotypes due to non-random linkage among alleles of loci not lying on the same chromosome. Two pseudo-haplotypes that imply a minor allele of TGF-βR2 or minor allele of VEGF-A SNPs combined with major alleles of the other SNPs were, respectively, associated with a protective effect, and susceptibility to NMSC. In addition, a pseudo-haplotype involving minor alleles of TGF-β2 rs900, TGF-βR1 rs334348 and rs4522809 SNPs might be a susceptibility marker for NMSC. In conclusion, our data suggest that a complex interplay among the genetic polymorphisms of TGF-β, TGF-β receptors and VEGF-A genes might influence the net effect of genetic background of the patients on NMSC development. This might be relevant in the risk evaluation, diagnosis and treatment of NMSC.


Introduction
The incidence of the vast majority of cancers, including melanoma, increases with advancing age. However, the mechanism by which aging influences the risk of developing cancer is not fully understood. Whether this is attributable to declines in immunity or other changes, such as changes in the aging tissue microenvironment, remains unclear [1].
Non-melanoma skin cancers (NMSC) include basal cell carcinoma (BCC) and squamous cell carcinoma (SCC), which are responsible for approximately 80 and 20% of all NMSC cases, respectively [2,3]. In humans, skin cancer arises in the epidermis or hair bulbs where proliferating basaloid tumor cells are detected. The extension and severity of these Patients were affected by non-melanoma skin cancer, clinically and histologically classified as basal cell carcinomas or squamous cell carcinomas. To perform comparative genotype analyses 161 age related healthy controls were also enrolled. Demographic and clinical characteristics of patients and control subjects are summarized in Table 1. In addition, both patients and controls belonged to same ethnic group, since their parents and grandparents were born in Western Sicily. Our study was performed in accordance with ethical standards of the Helsinki Declaration of the World Medical Association and Italian legislation, and was approved by the local institutional review board (Comitato Etico Palermo 1, protocol code CET1 03/2017, date of approval 1 March 2017). All participants gave their informed consent. Data were encoded to ensure privacy protection of patients and controls. Laboratory procedures were performed by lab staff unaware of the patient or control sources of the biological specimens.

Molecular Typing
As reported in Table 2, we selected six functional and common SNPs of the TGF-β pathway and one of the VEGF-A gene. Selection of these SNPs was performed using information acquired from dbSNP NCBI, the ENSEMBL database (http://www.ensembl. org/index.html, accessed on 3 April 2022), and the UCSC Genome Browser website (http: //genome.ucsc.edu, accessed on 3 April 2022). Alleles and genotypes were analyzed using on demand assays developed by KBioscience Ltd. (Middlesex, UK). Tests apply homogeneous Fluorescence Resonance Energy Transfer (FRET) detection and allele specific PCR (Kaspar) as previously described [30]. The genotypes were determined using the 7300 system SDS software, vs. 1.3 (Applied Biosystems, Monza (Mi) Italy) sample by sample, on the basis of the detection of unique (homozygous samples) or double (heterozygous samples) fluorescence signals.

Statistical Analysis
Allele and genotype frequencies were evaluated by gene count using an online statistical analysis tool applied to the evaluation of SNPs (https://www.snpstats.net/start.htm, accessed on 9 May 2022). Hardy-Weinberg equilibrium was checked by Pearson's test. Evaluation of power of calculation was performed with an online tool (https://sample-size.net/, accessed on 19 April 2022). Significant differences in genotype distributions between groups were calculated using chi-square or Fisher's exact test applied to the appropriate 2 × 2 or 3 × 2 contingency tables. Overdominant, codominant, dominant and recessive models of heredity were applied in multiple logistic regression analyses adjusting the results for age and gender. GraphPadInStat software (GraphPad, San Diego, CA, USA, version 3.06) and the above mentioned online statistical analysis tool applied to the evaluation of the association of SNPs with diseases in the presence of other variances (biological, genetic, or clinical), were used to calculate statistical significance (p-value cutoff < 0.05) of odds ratio (OR) values and 95% confidence intervals (95% CI). Epistatic interaction among the different SNPs was performed using the Haplotype frequencies estimation (frequency threshold for rare haplotype: 0.01) and Haplotype association with response estimation functions at https://www.snpstats.net/start.htm (accessed on 19 April 2022) searching for so called pseudo-haplotypes. Actually, the concept of epistasis analysis is very similar in theory to haplotype analysis [33,34], in particular when functionally well-defined pathways (i.e., TGF-β and VEGF functionally interacting pathways) are studied.

Allele and Genotype Frequencies of TGF-β and VEGF-A SNP in Non Melanocitic Skin Cancer Patient (NMSC) and Control Subjects
No significant differences were observed analyzing allelic and genotypic frequencies of TGF-β1, TGF-β2, and TGF-βR1 gene SNPs between NMSC patient and control subjects. When TGF-βR2 rs1800629 SNP was analyzed, significantly decreased frequencies of G allele and related genotypes were observed in the patient group (crude analysis) ( Table 3) These results were confirmed applying the dominant model of heredity in multiple logistic regression analyses, adjusting the results for age and gender (Table 3). In addition, both VEGF-A rs3025039T allele and TT homozygous genotype were found to be significantly increased in patients compared to the controls both in the crude analyses and when the recessive model of heredity in multiple logistic regression analyses was applied, adjusting the results for age and gender (Table 3).

SNP Frequencies of NMSC Patient and Control Groups Stratified According to 64 Years Age Cut Off and Gender
As reported in Table 1, the median age of the NMSC patients was 64 years old. This median value was chosen as cut-off for breakdown of the two patient and control populations in subgroups of <64 and ≥64 year old subjects. Data reported in Table 4 showed that by analyzing SNP frequencies stratified according to age cut off and adjusted by gender, there was an increase in TGF-βR2 rs4522809A homozygous genotype frequency in patients compared to the controls both in the <64 and ≥64 year old classes whereas heterozygous genotype and frequency of G positive genotypes on the whole were decreased.  On the other hand, when the analysis was applied to VEGF-A rs3025039 genotypes, significant differences were observed only between ≥64 years old patients and controls. In particular VEGF-A rs3025039CC homozygous genotype frequency was significantly decreased with respect to the control group, whereas heterozygous genotype and frequency of T positive genotypes on the whole was significantly increased. In addition, homozygous T genotype reached roughly a statistical significance when considered separately.
Data reported in Table 5 show that TGF-βR2 rs4522809A homozygous genotype frequency was significantly increased in patients respect to the controls both in females and in males. In addition, the frequency of G positive genotypes, on the whole, was significantly decreased (homozygous + heterozygous genotypes). Decreases of both homozygous and heterozygous G genotypes, considered separately, reached roughly statistical significance. On the other hand, when the analysis was applied to VEGF-A rs3025039 genotypes no significant differences were observed stratifying data according to gender, even if the frequency of rs3025039T homozygous genotype was increased with a statistical significancy near to p = 0.05 threshold (0.083 and 0.068, respectively, for female and male groups). The complete analyses of data stratified according to age cut-off and gender for all SNP types are reported in Supplementary Tables S1 and S2.

SNP Frequencies in NMSC Patients Stratified according to Skin Cancer Hystotype
SNP frequencies were analyzed stratifying patient data according to NMSC hystotype as affected by basal cell carcinoma (BCC) and as affected by squamous cell carcinoma (SCC) (adjusted by 64 years age cut off and gender), and compared with each other and with the control population (Table 6). This would mean that 20 patients were represented in the SCC and 50 in the BCC group. However, the analyses confirmed that TGF-β R2 rs4522809G positive genotype frequencies are reduced in NMSC patients with respect to the controls (BCC dominant model OR 0.24, CI 95% 0.12-0. 49

TGF-β and VEGF-A SNP Pseudo-Haplotype Associations with Non-Melanocytic Skin Cancer
Algorithms for haplotype frequencies estimation and association with the disease were used to evaluate the epistatic effect of combinations of genotypes at multiple loci on a complex trait such as NMSC due to the so-called gametic phase disequilibrium, as the alleles non-casually segregate within gametes, but are not physically located on the same chromosome [35]. This approach allows us to identify 31 functional pseudo-haplotypes composed by functional disequilibrium among alleles of functionally linked loci with a frequency > 0.01 in at least one of the populations analyzed, i.e., the NMSC or control population (Supplementary Table S4).
As reported in Table 7, pseudo-haplotype 2 composed by major alleles of TGF-β1, TGF-β2, TGF-βR1 and VEGF-A SNPs and minor allele of TGF-βR2 SNP showed a significantly reduced frequency in NMSC with respect to the control population. Table 7. Significantly different results obtained analyzing frequencies of 31 pseudo-haplotype (p-Hp) containing major or minor alleles of TGF-β and VEGF-A SNPs, in non-melanocytic skin cancer (NMSC) patients and control (CTRL) subjects. On the contrary, pseudo-haplotype 14, composed by major alleles of TGF-β1, TGF-β2, TGF-βR1 and TGF-βR2 SNPs and minor allele of VEGF-A had a significantly increased frequency in NMSC patients. In addition, pseudo-haplotype 24 composed by minor alleles of TGFβ2, and both TGFβR1 SNPs and major alleles of TGF-β1, TGF-β2 and VEGF-A SNPs showed a significantly increased frequency in NMSC respect to control subjects. Complete pseudo-haplotype statistical analyses are reported in Supplementary Table S4.

p-Hp
Taken together, these data seem to indicate that TGF-β and VEGF-A pathway modulation by gene variants might have a role in susceptibility to or protection against NMSC.

Discussion
Non-melanoma skin cancers, basal cell and squamous cell carcinoma (BCC and SCC) commonly arise in the epidermis or hair follicles where proliferating basaloid tumor cells are present. Cytokine and growth factors are involved in this process [36].
Vascular endothelial growth factors (VEGFs) are the key regulators in angiogenesis and have been shown to play a significant role in the progression and prognosis of angiogenesisrelated diseases, such as cancer [37][38][39][40]. VEGF is known to play a critical role in the development of non-melanoma skin cancers. Actually, VEGFs can affect skin carcinogenesis by directly interacting with keratinocytes, tumor cells, and immune cells [41,42]. In addition, VEGFs carry out other functions beyond their well-established effects on angiogenesis and influence skin carcinogenesis by directly activating tumor cells [19]. Moreover, several studies have uncovered direct effects of VEGF on keratinocytes and skin tumor cells. These studies suggested that in addition to enhancing angiogenesis, VEGF may promote skin carcinogenesis by altering the survival, proliferation, or stemness of keratinocytes and tumor cells in an autocrine manner [37,39,41,43].
Our data indicate that the frequency of VEGF-A rs3025039T positive genotypes was increased in NMSC patients, adjusting data for gender and age. Stratifying data according to age, T positive on the whole was increased in >64 years old subjects. Similarly, stratifying data according to gender, homozygous T genotype reached a roughly statistical significance when considered separately although this finding is dependent on a small number of subjects categorized according to gender and single genotype.
As above reported, the VEGF-A rs3025039T allele has been observed to give rise to lower circulating VEGF levels. However, studies on a South Italy population have demonstrated that the rs3025039C/T associated VEGF levels are variable and influenced by the presence of other VEGF-A SNP alleles in linkage disequilibrium with rs3025039C/T alleles [44]. A possible link between the rs3025039C/T SNP and several cancer types has been reported, and the contribution of this polymorphism to oncogenesis has been investigated in several types of tumors [45][46][47][48]. In some cases, T allele was found associated to protection against cancer [27,49]. Allele T associated susceptibility has been reported for oral cancer [50], head and neck cancer [51], breast cancer [52] and colorectal cancer [53,54]. In a recent paper, Kämmerer et al. [55], reporting the association of T allele with oral squamous cell carcinoma, showed that VEGF-A rs3025039T allele seems to be associated with an increased CD31 immuno-stain, currently used for measurement of microvessel density (MVD) and endothelial cell kinetic in normal and pathological tissues.
On the other hand, our data indicated that the statistical significance of the increase in the T positive genotype was maintained both in patients ≥64 years old and in patients affected by BCC. Recent findings suggest that in aging the reduced angiogenesis potential correlates with a reduced expression of VEGF-A [56]. Moreover, it appears intriguing that tumor cells of human BCCs tend to show weak VEGF expression with positive tumor cells predominantly localized to the invading margin [26,57]. In contrast, SCCs, which are typically more aggressive than BCCs, display more intense and widespread staining, with higher expression in tumor cells localized near infiltrating inflammatory cells [26]. However, our group of SCC patients is small as only 20 patients were represented in the SCC group, so the data obtained by comparing SCC patients with BCC and controls should be considered suggestive but not conclusive. In addition, the case mixes in both BCC and SCC groups (Table 1) do not allow us to further analyze possible significant differences associated with histotype, therapeutical approach, or relapses in relationship with SNP types.
Overall, the increased frequency of VEGF-A rs3025039T positive genotypes in our group of patients confirms the view that analysis of VEGFs genetic polymorphisms in aged cancer patients might be important for a better tailoring of cancer therapy [56].
The TGF-β signaling pathway plays an important role in the regulation of normal and cancer skin cell proliferation [13]. Increased expression of TGF-β and activation of the T regulatory pathway characterize the peri-tumoral skin microenvironment both in non-melanocytic tumors and in melanoma. TGF-β1 can modulate normal and pathological cell differentiation and proliferation in an auto-or paracrine manner [58]. In particular, the ability of TGF-β signaling to activate target genes enables the pathway to impact diverse cellular processes including proliferation, differentiation, migration, apoptosis, and endothelial cell remodeling [59,60] with different and sometimes opposite effects on cancer cells. Actually, TGF-β normally exerts tumor-suppressive effects. However, as the tumor advances in malignant progression, mutations in the TGF-β receptor system render the cancer cells unresponsive to the cytokine and to TGF-β through its multiple effects on non-malignant stromal cells, such as evasion of immune surveillance, creating a tumor-promoting factor that leads to increased invasion and metastasis [11].
In this scenario, given the role of genetic polymorphisms located in both transcribed and untranscribed regions of TGF-β and TGF-β receptors, the tuning of TGF-β effects might be crucial.
Accordingly, we found that TGF-βR2 SNP genotypes are differently distributed in NMSC patients and controls. In particular, TGF-βR2 rs4522809G positive genotype frequencies are reduced in NMSC patients with respect to the controls, adjusting data for gender and age. Similar results were obtained by analyzing genetic frequencies in female, male, younger or older patients compared to adequate controls. Moreover, this reduced frequency was observed both in SCC and BCC histologic types. Decreases of both homozygous and heterozygous G genotype reach roughly statistical significance when considered separately but this finding is dependent on a relatively small number of subjects categorized according to gender and single genotype.
It is not yet clear how TGF-β signaling is influenced by rs4522809. It has been hypothesized that the G allele of TGF-βR2 rs4522809, that is located in the intronic region, might modify the gene expression. Alternatively, it is possible that this gene variant regulates the amount and direction of transcription [61]. In this view, the increased frequency of rs4522809G positive genotypes in the control group allows to speculate about a possible protective effect of this gene variant against NMSC. It appears intriguing that a similar protective effect of the G allele was found in other types of cancer such as breast cancer [17]. We have also found that TGF-βR1 rs334349AA genotype frequencies are significantly increased in SCC patients compared both to BCC and controls. Association of this TGF-βR1 variant with cancer was documented in particular for colon rectal cancer [62]. Overall, our results, even if they need to be confirmed in a larger sample of NMSC patients (in particular for data comparing SCC and BCC histologic types) suggest that the net effect of the TGF-β gene pathway on skin cancer pathogenesis is influenced by polymorphisms of genes coding for the TGF-β receptor complex.
As is well known, the net effect of the presence of a functionally relevant SNP might be masked or counterbalanced by other SNPs in the same or in interacting genes [31]. Our data in their complexity seem to suggest both protective and susceptible polymorphism in strictly interacting genes as TGF-β family and VEGF-A [63,64]. In this regard, we have explored the interaction among the six SNP types, applying an online haplotype searching tool to determine the so called functional pseudo-haplotypes [32]. This would describe the correlation among genotypes of loci not lying on the same chromosome, that has been also described as "gametic phase disequilibrium", i.e., the non-random association of alleles within gametes, even for physically unlinked loci tethered on different chromosomes [65]. This approach allowed us to establish that two pseudo-haplotypes showed a significantly different distribution frequency in NMSC with respect to the control population, implicating a minor allele of TGF-βR2 or minor allele of VEGF-A SNPs combined with major alleles of the other SNPs studied. Analyzing the frequencies of the two TGF-βR2 and VEGF-A SNPs individually, it seems that the first haplotype is associated with a protective effect against NMSC and the second with susceptibility. These findings would suggest that the effect of TGF-βR2 and VEGF-A SNPs might be more evident when the major alleles of the other loci are present. It is intriguing to observe that when minor alleles of TGF-βR2 and VEGF-A SNPs are contemporaneously associated in the same pseudo-haplotype, no significant differences between frequencies of NMSC and control were observed (Table  S4 pseudo-haplotypes 12 and 16). Finally, pseudo-haplotype analyses revealed that the pseudo-haplotype 24 composed by minor alleles of TGF-β2 rs900, TGF-βR1 rs334348 and rs4522809 SNPs and major alleles of TGF-β1, TGF-βR2 and VEGF-A SNPs might be a susceptibility factor for NMSC.
On the other hand, both next generation sequencing studies and GWAS have not directly evidenced that SNP of TGF-β pathway or VEGF-A are genetic markers of risk for NMSC. These approaches have highlighted that pigmentation genes as TYR, TYRP1, OCA2, SLC24A5, SLC45A2, POMC, ASIP, ATRN, involved in activation of the Hedgehog (HH) pathway, were identified as susceptible loci for BCC and SCC, independently of the pigmentation phenotypes [66,67]. HLA class II, TP63, FOXP1 genes were associated with NMSC [68,69]. However, novel approaches that integrate skin expression-related singlenucleotide polymorphisms (eSNPs) and pathway analysis based on Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), or BioCarta databases allowed identification of the role of genes of the TGF-β pathway interacting with genetic activation of the Hedgehog pathway [70,71]. The NO pathway in the BioCarta includes VEGFA SNPs and plays an important role in the regulation of vascular endothelial function facilitating, in particular, BCC development and an angiogenic response propagated by the growing skin cancer [72]. In this view, it is hypothesized that VEGF-A and TGF-βR2 genetic variants studied by our group can be considered components of multiloci and complex metabolic pathways that characterize NMSC development.

Conclusions
In conclusion, our data suggest that a complex interplay among the genetic polymorphisms of TGF-β, TGF-β receptors and VEGF-A genes might influence the net effect of the genetic background of the patients on development of non-melanocytic tumors of the skin. Actually, in humans, their gene products can mediate stimulating or inhibiting cell growth effects depending on the cellular and tissue targets [10]. TGF-β, in particular, might mediate different and time-differentiated effects on tumor evolution. The apparent contrasting effects might be due to progressively distinct genetic states and modifications of the tumor microenvironment influencing proliferation, metabolism and expression of key genes targeted [34,58,73]. In this view, larger investigations are warranted to explore the relevance of the TGF-β and VEGF pathway genetic background as a diagnostic or predictive marker that can be targeted to treat non-melanocytic skin cancers.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13071235/s1, Table S1: SNP frequencies in Non Melanocitic Skin Cancer Patients (NMSC) and Controls (CTRL) stratified according to 64 years age cut off (adjusted by gender); Table S2: SNP frequencies in Non Melanocitic Skin Cancer Patients (NMSC) and Controls (CTRL) stratified according to gender (adjusted by age cut off); Table S3: SNP frequencies in Non Melanocitic Skin Cancer Patients (NMSC) affected by basal cell carcinoma (BCC) compared to NMSC with squamous cell carcinoma (SCC) (adjusted by 64 years age cut off and gender); Table S4: Analysis of frequencies of 31 pseudo-haplotype (p-Hp) containing major or minor alleles of TGF-β gene family and VEGF-A SNPs, in non-melanocytic skin cancer (NMSC) patients and controls (CTRL).  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: All data generated or analyzed during this study are stored in electronic archives and can be supplied on request.