The Contribution of Genetic Variants to the Risk of Papillary Thyroid Carcinoma in the Kazakh Population: Study of Common Single Nucleotide Polymorphisms and Their Clinicopathological Correlations

Objective: Risk for developing papillary thyroid carcinoma (PTC), the most common endocrine malignancy, is thought to be mediated by lifestyle, environmental exposures and genetic factors. Recent progress in the genome-wide association studies of thyroid cancer leads to the identification of several genetic variants conferring risk to this malignancy across different ethnicities. We set out to elucidate the impact of selected single nucleotide polymorphisms (SNPs) on PTC risk and to evaluate clinicopathological correlations of these genetic variants in the Kazakh population for the first time.

Methods: Eight SNPs were genotyped in 485 patients with PTC and 1,008 healthy control Kazakh subjects. The association analysis and multivariable modeling of PTC risk by the genetic factors, supplemented with rigorous statistical validation, were performed.

Result: Five of the eight SNPs: rs965513 (FOXE1/PTCSC2, P = 1.3E-16), rs1867277 (FOXE1 5’UTR, P = 7.5E-06), rs2439302 (NRG1 intron 1, P = 4.0E-05), rs944289 (PTCSC3/NKX2-1, P = 4.5E-06) and rs10136427 (BATF upstream, P = 9.8E-03) were significantly associated with PTC. rs966423 (DIRC3, P = 0.07) showed a suggestive association. rs7267944 (DHX35) was associated with PTC risk in males (P = 0.02), rs1867277 (FOXE1) conferred the higher risk in subjects older than 55 years (P = 7.0E-05), and rs6983267 (POU5F1B/CCAT2) was associated with pT3–T4 tumors (P = 0.01). The contribution of genetic component (unidirectional independent effects of rs965513, rs944289, rs2439302 and rs10136427 adjusted for age and sex) to PTC risk in the analyzed series was estimated to be 30–40%.

Conclusion: Genetic factors analyzed in the present work display significant association signals with PTC either on the whole group analysis or in particular clinicopathological groups and account for about one-third of the risk for PTC in the Kazakh population.

Introduction

Papillary thyroid carcinoma (PTC), a well-differentiated thyroid cancer of follicular cell origin, accounts for about 80% of all thyroid cancers worldwide being the most common endocrine malignancy. According to the IARC, in 2018 thyroid cancer affected 567,233 patients worldwide, making it the 9^th most prevalent human cancer (3.1%) with the average age-standardized incidence of 6.7 and mortality rate of 0.42 per 100,000 of population (1). Region-specific incidence rates vary broadly from 1.0 in Micronesia to 15.0 in North America per 100,000 of population. In Kazakhstan, the age-standardized incidence of thyroid cancer was 2.4 per 100,000 of population accounting for 1.4% of all newly diagnosed cancers in the country in 2018.

With the improvements in cancer detection and diagnosis, the incidence of thyroid cancer is growing in most countries displaying one of the fastest increases in rate among common cancers. While the advances in thyroid nodule visualization such as ultrasound imaging and their facile assessment using ultrasound-guided fine-needle aspiration cytology have likely contributed to this uptrend, the additional reasons for the increasing incidence are investigated. Besides of well-established risk factor for thyroid cancer such as ionizing radiation, other environmental agents, including iodine deficiency, natural and technogenic pollutants with hormone disrupter effects, exposures to excessive nitrate (2, 3) and various chemicals are discussed or considered.

As a complex disease, PTC is thought to be dependent not only on environmental, but also on genetic factors. Studies of familial thyroid cancer estimated the contribution of genetic component to the risk of disease to be ranging from 28 to 53% (4, 5). At the population level, hereditary factors possibly contributing to the phenotype (e.g. the development of a condition or a disease) are usually identified in genetic association studies. To date, a number of well-powered genome-wide association studies (GWAS) or target gene investigations in thyroid cancer have been performed in the groups of different ethnicities in non-exposed or exposed to radiation individuals (6–14). GWAS findings and consequent independent replication studies have convincingly demonstrated robust associations of rs965513 (FOXE1, forkhead box E1 and/or PTCSC2, papillary thyroid carcinoma susceptibility candidate 2; chromosome 9q22.33), rs944289 (PTCSC3, papillary thyroid carcinoma susceptibility candidate 3 and/or NKX2-1, NK2 homeobox 1; 14q13.3), rs1867277 (FOXE1; 9q22.33), rs2439302 (NRG1, neuregulin 1; 8q12) and rs966423 (DIRC3, disrupted in renal carcinoma 3; 2q35) SNPs with differentiated thyroid cancer, principally PTC (15–27), reviewed in (28). The strength of association signal for these SNPs in terms of odds ratios (OR) ranged from 1.28 to 1.70 in most studies. More recent studies have identified associations between the rs6983267 (POU5F1B, POU class 5 homeobox 1B and/or CCAT2, colon cancer associated transcript 2; 8q24) and thyroid cancer in different populations. A systematic review with meta-analysis of four studies that included a total of 2,825 cases and 9,684 controls confirmed the G allele of the rs6983267 to be a risk factor for thyroid cancer with an OR = 1.08, P = 0.01 (29).

The novel GWAS candidate loci continue to emerge. A recent combined analysis of GWAS results and the Italian replication study provided evidence of association of risk for differentiated thyroid cancer with rs10136427 (BATF, basic leucine zipper ATF-like transcription factor, 14q24.3) with an OR = 1.40, P = 4.35E-07) and rs7267944 (DHX35, DEAH-box helicase 35, 20q12) with an OR 1.39, P = 2.13E-08. These associations were replicated in the Polish and Spanish populations with little evidence of population heterogeneity (the combined, OR = 1.30, P = 9.30E-07 and OR 1.32, P=1.34E-08, respectively) (10).

To the best of the authors’ knowledge, studies of rs10136427 (BATF, 14q24.3) and rs7267944 (DHX35, 20q12) in PTC have not been replicated in independent studies. We therefore aimed to examine the six well-described SNPs discussed above, namely rs965513 (FOXE1/PTCSC2, 9q22.33), rs944289 (PTCSC3/NKX2-1, 14q13.3), rs1867277 (FOXE1; 9q22.33), rs2439302 (NRG1, 8q12), rs966423 (DIRC3, 2q35) and rs6983267 (POU5F1B/CCAT2, 8q24.2), and two SNPs newly discovered to be associated with thyroid cancer (10), rs10136427 (BATF, 14q24.3) and rs7267944 (DHX35, 20q12) in a relatively large case–control series. This work is the first to characterize the eight SNPs in the Kazakh population. In addition to the classical association analysis, we estimated the contribution of the genetic variants to PTC risk, and assessed the relationships with clinicopathological characteristics of tumors in the study since available information is very limited in the literature.

Material and Methods

Study Population

A total of 485 patients with histologically confirmed PTC (90.3% females, mean age 54.78 ± 13.3 y.o., 18–87 y.o., range) operated from 1980 to 2015, and 1,008 healthy control subjects (78.7% females, mean age 39.0 ± 15.8 y.o., 15–83 y.o., range) of Kazakh origin were recruited. Clinicopathological information was retrieved from medical records (Table 1). The pathological classification was based on the World Health Organization definitions (30), pathological staging (pTNM, where the T category defined the anatomic extent of cancer for the tumor, N for the lymph nodes and M for distant metastases) was according to the 7th edition of TNM classification system (31). Patients and control subjects had no history of radiation exposure. All participants or their parents/guardian gave written informed consent in accordance with the Declaration of Helsinki. A peripheral venous blood sample was collected from each participant. The protocol of this study was approved by corresponding ethics committees.

TABLE 1

Table 1 Demographic characteristics of control subjects and PTC patients, and clinicopathological data.

DNA Isolation and Genotyping

Blood DNA was extracted using QIAamp DNA Mini Kit (QIAGEN, Tokyo, Japan) according to the manufacturer’s protocol. DNAs of sufficient quality for genotyping were obtained from all 485 PTC patients and 1,008 control participants.

Genotyping was performed with predesigned Custom Applied Biosystems TaqMan SNP Genotyping Assays (Table 2) using TaqMan Genotyping Master mix (all reagents from ThermoFisher Scientific) and 10 ng genomic DNA per 10 µl reaction in a Light Cycler 480 (Roche, Indianapolis, IN). Cycling conditions were as follows: denaturation at 95°C for 10 min followed by 60 cycles of 92°C for 15 s and 62°C for 1 min for all SNPs. As a quality control, 15–20% of all samples were randomly selected and re-run in duplicates for each SNP. Full concordance between the experiments was observed.

TABLE 2

Table 2 TaqMan primer/probe set used for genotyping.

Statistical Analyses

Association Analysis

We used PLINK 1.9 (32) software to run the multiplicative genetic models in the case–control sample for each SNP with age and sex as covariates. This type of model evaluates the impact of individual alleles of a polymorphic locus on the disease. The multiplicative models have been used in the vast majority of the genome-wide and replication association studies of thyroid cancer (6–27); using those in our work provided an opportunity to compare the strength of association signals (ORs) between the previous studies and our findings. The risk alleles were assigned according to the cited literature sources for consistency; summary information on the risk alleles is provided in Table 3. Multiple testing correction (the Benjamini–Hochberg method) and the adaptive label-swapping permutation test (10⁶, maximum) were performed using options available in PLINK.

TABLE 3

Table 3 Summary information on the risk alleles of analyzed SNPs.

Associations between each SNP and clinicopathological parameters of PTCs were assessed using logistic regression models with binary outcomes sex (F vs M), age (≥55 vs <55 years old), pathological tumor (pT) category (pT3 or pT4 vs pT1 or pT2) or nodal disease (N1 vs N0, i.e. present vs absent) as dependent variables, and individual SNPs, age and sex (where applicable) as explanatory variables. The LOGISTIC procedure in the 3.71 release of SAS Studio for the 9.4M5 version of SAS (SAS Institute, Cary, NC, USA) was used for these calculations.

Exact two-sided tests, permutation tests and exact test for equality of allele frequencies for stratified groups were performed using the ‘HardyWeinberg’ package in R (37).

All tests were two-sided, P <0.05 was considered statistically significant.

Predictive Modeling of Papillary Thyroid Carcinoma

To evaluate the impact of the genetic component on PTC risk in the given case–control sample, we used multivariable logistic regression modeling. The initial full model included all eight SNPs in the study, and age and sex as explanatory variables. The reduced model was determined by stepwise or non-automatic selection of variables to achieve minimum Akaike information criterion. Statistical validation of the reduced model was performed using permutation analysis as described before (38), and bootstrapping with 0.9 sampling rate (i.e., selecting 90% of data for each sample using the unrestricted random sampling method) in 10,000 replicates using the SURVEYSELECT procedure. The receiver operating characteristic (ROC) analysis was performed to assess the predictive performance of the reduced model, supplemented with the leave-one-out cross-validation.

Results

Single_Nucleotide Polymorphisms Association With and Impact on Papillary Thyroid Carcinoma Risk

Five of the eight SNPs displayed significant association signals in the Kazakh population with ORs similar to those in the original studies and follow-up publications (Table 4). The strongest associations between a risk allele and sporadic PTC were observed for rs965513 (FOXE1/PTCSC2, 9q22.33; OR = 2.25, P = 1.3E-16), rs1867277 (FOXE1 5’UTR, 9q22.33; OR = 1.52, P = 7.5E-06), rs2439302 (NRG1 intron 1, 8q12; OR = 1.46, P = 4.0E-05), rs944289 (PTCSC3/NKX2-1, 14q13.3; OR = 1.44, P = 4.5E-06), and rs10136427 (intergenic region upstream BATF, 14q24.3; OR = 1.30, P = 9.8E-03). rs966423 (DIRC3, 2q35) showed a significant association (OR = 1.25, P = 5.8E-03) on unadjusted analysis, but significance became marginal after adjusting for age and sex (OR = 1.18, P = 0.07). Adjustment for multiple testing and statistical validation (permutation) confirmed significant association of the five SNPs (rs965513, rs944289, rs1867277, rs2439302, rs10136427) and suggestive association for rs966423.

TABLE 4

Table 4 Association analysis of papillary thyroid carcinoma in the Kazakh population.¹

Two remaining SNPs, rs7267944 (DHX35, 20q12; OR = 1.04, P = 0.71) and rs6983267 (POU5F1B/CCAT2, 8q24.2; OR = 1.09, P = 0.36) did not display significant associations at this stage of analysis.

After obtaining evidence that certain examined SNPs display statistically significant association signals, we set out to determine the performance of a statistical model of the risk for PTC based on genetic factors. The reduced model included four SNPs: rs965513 (FOXE1/PTCSC2), rs944289 (NKX2-1, PTCSC3), rs2439302 (NRG1) and rs10136427 (BATF) (Table 5). Their association signals remained significant after correction for multiple testing (Bonferroni and FDR). Statistical validation confirmed significant association of these SNPs with cancer risk (permutation), and confidence intervals almost did not change on bootstrapping. Of note, the OR estimates for the four SNPs in the model were very similar to those obtained in the single-SNP models (see Table 4) indicative of independent contribution of each SNP to thyroid cancer risk. The area under the ROC curve (AUC) was 0.82 (95%CI 0.80–0.84; P = 3.2E-183 as compared with a random classifier); cross-validation did not demonstrate model overfit returning a similar AUC of 0.82 (0.80–0.84). The Cox and Snell pseudo-R² of the optimal model was 0.27, and the Nagelkerke pseudo-R² was 0.38 suggesting that some 30–40% of variability in the risk for PTC in the analyzed series could be explained by the model that included age, sex and the four SNPs.

TABLE 5

Table 5 The SNP-based logistic regression model of risk for PTC in the Kazakh population.¹

Single-Nucleotide Polymorphism Association With Clinicopathological Parameters of Papillary Thyroid Carcinoma

We assessed the genetic variants association with patients’ sex, age (older than 55 years old), the pT category (pT3–T4 vs pT1–T2), and regional metastasis (N1). As shown in Table 6, very few statistically significant associations were found. rs7267944 (DHX35) was negatively associated with female sex (OR = 0.40, P = 6.1E-05), rs1867277 (FOXE1) was weakly but nominally significantly associated with the older patients’ age (OR = 1.32, P = 0.03) and rs6983267 (POU5F1B/CCAT2) with more advanced tumors (OR = 1.49, P = 0.01).

TABLE 6

Table 6 Association of the genetic variants with clinicopathological characteristics.¹

The strong association of rs7267944 (DHX35) with patients’ sex prompted us to test its association with PTC risk using stratified sampling. While no association was found in females OR = 0.95 (95%CI 0.76–1.18, P = 0.62 adjusted for age), the association signal in males was significant with an OR = 1.83 (95%CI 1.09–3.09, P = 0.02 adjusted for age). The difference in effect size was statistically significant (P = 0.023, the Breslow–Day test). No deviations from Hardy–Weinberg equilibrium in the groups of PTC patients or healthy control subjects, either non-stratified or stratified by sex, was found for this genetic variant (P > 0.4 for any exact two-sided test, P > 0.4 for any permutation test). Exact test for equality of allele frequencies for males and females in the control subjects was negative (P = 0.06), but in PTC patients a strong inequality was observed (P = 8.26E-05), in line with other statistical findings.

The modifying effect of age on rs1867277 (FOXE1) was tested in respective groups of patients and control subjects younger or older 55 years old. In the younger group, rs1867277 displayed an association signal with OR = 1.44 (95%CI 1.16–1.80, P = 9.5E-04 adjusted for sex), and in the older group with OR = 1.84 (95%CI 1.36–2.49, P = 7.0E-05 adjusted for sex). There was no surprise that the association was significant in both age groups as rs1867277 was significant on the whole group association analysis, which was adjusted for age (see Table 4). Clearly, the effect of rs1867277 on the risk for PTC was more pronounced in subject older than 55 years old, although the difference did not reach statistical significance (P = 0.20). No deviations from Hardy–Weinberg equilibrium in the groups of PTC patients or healthy control subjects were found for rs1867277 (P > 0.1 for any exact two-sided test, P > 0.1 for any permutation test). Exact test for equality of allele frequencies for subjects younger than 55 years old was negative (P = 0.85), while in older subjects the inequality existed (P = 0.03).

We also tested the association of rs6983267 (POU5F1B/CCAT2) with PTC of different pT stage. In pT1–T2 tumors, the association was insignificant with OR = 1.0 (95%CI 0.82–1.21, P = 0.96 adjusted for age and sex), while in pT3–T4 tumors the signal was significant, OR = 1.47 (95%CI 1.09–1.98, P = 0.01 adjusted for age and sex). The difference in ORs was statistically significant (P = 0.03). No deviations from Hardy–Weinberg equilibrium in the groups of PTC patients or healthy control subjects were found for rs6983267 (P > 0.1 for any exact two-sided test, P > 0.1 for any permutation test). Exact test for equality of allele frequencies in the subgroup of PTCs of pT1–pT2 category was negative (P = 0.32) and marginally significant in pT3–pT4 PTCs (P = 0.07).

Discussion

In the present study we set out to determine the impact of genetic factors on PTC risk in the Kazakh population. We focused on several SNPs found to display confident association signals in the previous studies across different populations/ethnicities and also tested two SNPs that have been newly discovered in GWAS of thyroid cancer.

Our genotyping results unambiguously confirmed the associations of rs965513 (FOXE1/PTCSC2, 9q22.33), rs1867277 (FOXE1 5’UTR, 9q22.33), rs944289 (PTCSC3/NKX2-1, 14q13.3), rs2439302 (NRG1 intron 1, 8q12) using canonical multivariate analysis essentially supplemented by rigorous statistical validation.

Functional roles of rs965513 and rs1867277 located in the FOXE1 locus on chromosome 9q22.33 were linked to the transcriptional regulation of FOXE1 and PTCSC2. rs965513 was shown to affect the expression of FOXE1, PTCSC2 and TSHR (thyroid stimulating hormone receptor) in thyroid tissue (34), and rs1867277 regulates FOXE1 expression through the recruitment of USF1/USF2 transcription factors (15), implicating these SNPs into thyroid homeostasis and development. Of note, transgenic mice overexpressing FOXE1 in their thyroids displayed retardation in the proliferation of follicular cells, suggestive of its tumor suppressor function (39). A meta-analysis that combined data from 23 studies in different countries and ethnicities evaluated that rs965513[A] risk allele had an OR = 1.58 (95% CI 1.32–1.90) in the pooled sample, and OR = 1.65 and 1.49 in Caucasian and Asian populations, respectively (40). Interestingly, in the Kazakh population, which is of Asian descent, we found an OR = 2.25 (95% CI 1.85–2.73), which is one of the strongest association signals ever reported for the FOXE1 locus. Also of interest is the finding of age-related effect rs1867277 (i.e., the higher risk in patients aged more than 55 years), which is reported for the first time. It should be noted that despite rs965513 and rs1867277 are located in the same genetic locus, their effect on PTC risk is independent (41). The age relatedness of rs1867277 effect could be addressed in independent or already available studies.

rs944289 located on chromosome 14q13.3 regulates expression of the PTCSC3 lncRNA, which has tumor suppressor effect in thyroid cancer cell lines, through the recruitments of C/EBPα and β transcription factors (35). PTCSC3 level was found to be significantly lower in PTC than in normal thyroid tissue (26), which corresponds well with its tumor suppressor function. Interestingly, rs944289 besides of PTC was also associated with follicular adenoma (26), indicating its broader function in thyroid tumorigenesis.

rs2439302 is located in intron 1 of NRG1 on chromosome 8q12. NGR1 encodes human epidermal growth factor receptor 3 (HER3) ligand whose dimers with HER2 can activate MAPK and AKT pathways known to play an important role in PTC (42). Similarly to PTCSC3, NRG1 was also earlier associated with follicular adenoma (26).

rs10136427, whose association with PTC risk was confirmed in the Kazakh population for the first time, was previously detected on GWAS in the Italian, Polish, and Spanish population study providing strong evidence of association with DTC (12). rs10136427 is located in an intergenic region upstream BATF. BATF proteins are the “AP-1 inhibitors”; findings in mouse myeloid leukemia cells suggested they can act as tumor suppressors by promoting cell growth arrest and cell differentiation. Whether BATF could play a similar role in other tissues, such as the thyroid, remains unknown. Since this genetic variant was not immediately associated with BATF expression, a possibility was suggested that this genetic locus may act as a trans-regulatory region controlling the expression of distant genes that reside on the same or even different chromosome(s) (trans-eQTL) (12).

rs966423 located in the DIRC3 intron on chromosome 2q35, displayed marginally significant association (P = 0.07) in the Kazakh population. This genetic variant was significantly associated with the risk for thyroid cancer in both European and Asian ethnicities with ORs from 1.28 to 1.34 (reviewed in (28). In our study the OR = 1.18 (adjusted for age and sex) is lower than those previously described. It therefore is likely that our sample size (485 PTC patients and 1,008 healthy control subjects) did not provide sufficient statistical power (achieved power 44%). We interpret the association signal of this genetic variant as suggestive in the Kazakh population. The functional role of DIRC3 lncRNA is likely to be that of tumor suppressor, and its relevance to thyroid cancer (8, 9) and other human malignancies such as, originally, familial renal cancer (43), melanoma (44), breast cancer (45) and laryngeal squamous cell carcinoma (46) has been reported.

The finding for rs7267944, which is located approximately 280 kB telomeric to DHX35 on chromosome 20q12, was somewhat unexpected. While on the whole group association analysis the signal of rs7267944 was insignificant, we noticed a strong modifying effect of sex. Accordingly, we found significant association of rs7267944 with PTC in males but not in females. DHX35 encodes a putative RNA helicase of DEAD/DEAH-box family, which are implicated in translation initiation, RNA splicing, and ribosome and spliceosome assembly. DHX35 is relatively highly expressed in the endometrium, ovaries, prostate and testes possibly pointing at its relatedness to sex-specific biological function (47). In the thyroid, DHX35 is also expressed, although its role in tissue homeostasis and carcinogenesis remains unestablished. Within our study we could not determine the reason for DHX35 association with sex, which, besides the true association could be sampling bias, a phenomenon specific for the given ethnic group (and relevant environmental exposures) or occurring by chance. This could be clarified in an independent study in the Kazakh population and also in other ethnic groups by researcher with access to rs7267944 [or other SNP(s) in strong or perfect linkage disequilibrium with it] genotyping data and clinical/demographic information. It is also possible that rs7267944 may point on the genetic factor other than DHX35 on chromosome 20q12 or elsewhere due to trans-eQTL effect.

A recent GWAS has identified rs6983267 (POU5F1B/CCAT2) as a key locus in the 8q24 region previously associated with DTC/PTC. However, the association with thyroid cancer was somewhat controversial since the significant association was found in the Polish and UK populations, but no association was found in the Spanish, Italian, and Japanese groups (10). In the Kazakh population under study, we did not observe significant association signal on the whole group analysis, yet a correlation with the higher pT tumor stage was detected. When pT3–T4 tumors were tested, a significant association with PTC risk was confirmed. It is tempting to relegate controversies in the rs6983267 association with thyroid cancer in different populations not only to genetic heterogeneity but also to different distribution of clinicopathological characteristics of tumors in country-specific samples. rs6983267 resides in the intronic region of POU5F1B, which encodes a transcriptional activator implicated in multisite cancers (48–51). It is worth noting that this genetic variant is also localized inside the CCAT2 lncRNA upregulated in colon cancer and implicated in other human malignancies (52–55). The exact roles of either POU5F1B or CCAT2 in PTC remain to be established.

After confirming the associations of rs965513 (FOXE1/PTCSC2), rs944289 (PTCSC3/NKX2-1), rs1867277 (FOXE1), rs2439302 (NRG1), rs10136427 (BATF) and, suggestively, of rs966423 (DIRC3) with PTC in the whole group or on subgroup analysis for rs7267944 (DHX35) and rs6983267 (POU5F1B/CCAT2), we combined these genetic variants in a statistical model to evaluate their contribution to PTC risk as the predictive strength of the genetic variants can be improved by combining multiple SNPs in a model (27, 56). The final model, which included four SNPs, rs965513 (FOXE1/PTCSC2), rs944289 (NKX2-1, PTCSC3), rs2439302 (NRG1) and rs10136427 (BATF), was subjected to statistical validation to ensure its reliability. The model had good predictive strength as judged by the ROC analysis (AUC = 0.82). Using two different analogs of the coefficient of determination for logistic regression models, we considered it safe to claim that the four SNPs in the optimal model, adjusted for age and sex effects, could explain about 30–40% of the risk for PTC in the Kazakh population examined in the study with a retrospective case–control design.

Our study had certain advantages such as homogenous ethnicity of the participants, large sample size that provided sufficient statistical power to detect meaningful associations, and thorough selection of the genetic variants. On the other hand, the study was not devoid of limitations. We could not fully rule out sampling bias and acknowledge insufficient age and sex matching of cases and controls, which could affect the accuracy of some of the results obtained in the study. Also, clinicomorphological information was not detailed enough to analyze potentially clinically relevant correlations, and the lack of data on the participants’ lifestyles and environmental exposures impeded the assessment of the impact of these factors on PTC risk.

In summary, our results unambiguously demonstrate the existence of genetic determinants of susceptibility to PTC among the SNPs analyzed in this work in the Kazakh population. We confirm the associations of rs965513 (FOXE1/PTCSC2), rs944289 (PTCSC3/NKX2-1), rs1867277 (FOXE1), rs2439302 (NRG1), and rs10136427 (BATF). The association of rs966423 (DIRC3) with PTC risk in the Kazakh population is suggestive. The association signals in terms of ORs were generally comparable to those in typical Asian and European populations, and that of rs965513 (FOXE1/PTCSC2) was the highest so far reported. We also detected the age-related effect of rs1867277 (FOXE1) conferring the higher risk for PTC in patients older than 55 years and the association of rs7267944 (DHX35) with PTC risk in males and that of rs6983267 (POU5F1B/CCAT2) with more advanced tumor pT stage. We estimate the contribution of genetic factors to the susceptibility to PTC in the analyzed series from Kazakhstan to 30–40%, accounting for age and sex. To better understand the impact of different factors affecting PTC risk, further studies would be desirable to increase the number of potential genetic loci and to include the data on individual lifestyle and exposures to environmental agents.

Data Availability Statement

Dataset used in this study may be made available upon request to Dr. Zhanna Mussazhanova (mussazhanova@nagasaki-u.ac.jp), subject to approval by the Semey Medical University Ethics Committee. Information on all genetic variants analyzed in this study is available in NCBI dbSNP database (https://www.ncbi.nlm.nih.gov/snp/).

Ethics Statement

The studies involving human participants were reviewed and approved by the Semey Medical University Ethics Committee and Nagasaki University Human Genome Ethics Committee. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

ZM, VS, ME, ZKa, RI, SY, and MN conceived and designed the study. ZM, AK, ME, ZKa, RI, ZY, MM, AM, MS, ZKo, and SB collected and formalized the clinical information. ZM and TR performed the experiments, analyzed and formalized the raw data. VS and HK performed the statistical analyses. TR, SY, KM, and MN contributed comments to the paper. ZM and VS wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported in part by the Atomic Bomb Disease Institute of Nagasaki University, The Joint Hiroshima University, Nagasaki University, and Fukushima Medical University Research Base for Radiation Accidents and Medical Science, and by Takeda Science Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank all participants of this study for their collaboration.

References

1. International Agency for Research on Cancer. Cancer Today (powered by GLOBOCAN) (2019). Available at: https://gco.iarc.fr/today/home (Accessed March 12, 2020).

Google Scholar

2. Drozd VM, Saenko VA, Brenner AV, Drozdovitch V, Pashkevich VI, Kudelsky AV, et al. Major Factors Affecting Incidence of Childhood Thyroid Cancer in Belarus after the Chernobyl Accident: Do Nitrates in Drinking Water Play a Role? PLoS One (2015) 10:e0137226. doi: 10.1371/journal.pone.0137226

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ward MH, Kilfoy BA, Weyer PJ, Anderson KE, Folsom AR, Cerhan JR. Nitrate intake and the risk of thyroid cancer and thyroid disease. Epidemiology (2010) 21:389–95. doi: 10.1097/EDE.0b013e3181d6201d

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Czene K, Lichtenstein P, Hemminki K. Environmental and heritable causes of cancer among 9.6 million individuals in the Swedish Family-Cancer Database. Int J Cancer (2002) 99:260–6. doi: 10.1002/ijc.10332

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kerber RA, O’Brien E. A cohort study of cancer risk in relation to family histories of cancer in the Utah population database. Cancer (2005) 103:1906–15. doi: 10.1002/cncr.20989

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gudmundsson J, Sulem P, Gudbjartsson DF, Jonasson JG, Sigurdsson A, Bergthorsson JT, et al. Common variants on 9q22.33 and 14q13.3 predispose to thyroid cancer in European populations. Nat Genet (2009) 41:460–4. doi: 10.1038/ng.339

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Takahashi M, Saenko VA, Rogounovitch TI, Kawaguchi T, Drozd VM, Takigawa-Imamura H, et al. The FOXE1 locus is a major genetic determinant for radiation-related thyroid carcinoma in Chernobyl. Hum Mol Genet (2010) 19:2516–23. doi: 10.1093/hmg/ddq123

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Gudmundsson J, Sulem P, Gudbjartsson DF, Jonasson JG, Masson G, He H, et al. Discovery of common variants associated with low TSH levels and thyroid cancer risk. Nat Genet (2012) 44:319–22. doi: 10.1038/ng.1046

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kohler A, Chen B, Gemignani F, Elisei R, Romei C, Figlioli G, et al. Genome-wide association study on differentiated thyroid cancer. J Clin Endocrinol Metab (2013) 98:E1674–81. doi: 10.1210/jc.2013-1941

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Figlioli G, Kohler A, Chen B, Elisei R, Romei C, Cipollini M, et al. Novel genome-wide association study-based candidate loci for differentiated thyroid cancer risk. J Clin Endocrinol Metab (2014) 99:E2084–92. doi: 10.1210/jc.2014-1734

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Mancikova V, Cruz R, Inglada-Perez L, Fernandez-Rozadilla C, Landa I, Cameselle-Teijeiro J, et al. Thyroid cancer GWAS identifies 10q26.12 and 6q14.1 as novel susceptibility loci and reveals genetic heterogeneity among populations. Int J Cancer (2015) 137:1870–8. doi: 10.1002/ijc.29557

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Figlioli G, Chen B, Elisei R, Romei C, Campo C, Cipollini M, et al. Novel genetic variants in differentiated thyroid cancer and assessment of the cumulative risk. Sci Rep (2015) 5:8922. doi: 10.1038/srep08922

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gudmundsson J, Thorleifsson G, Sigurdsson JK, Stefansdottir L, Jonasson JG, Gudjonsson SA, et al. A genome-wide association study yields five novel thyroid cancer risk loci. Nat Commun (2017) 8:14517. doi: 10.1038/ncomms14517

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Son HY, Hwangbo Y, Yoo SK, Im SW, Yang SD, Kwak SJ, et al. Genome-wide association and expression quantitative trait loci studies identify multiple susceptibility loci for thyroid cancer. Nat Commun (2017) 8:15966. doi: 10.1038/ncomms15966

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Landa I, Ruiz-Llorente S, Montero-Conde C, Inglada-Perez L, Schiavi F, Leskela S, et al. The variant rs1867277 in FOXE1 gene confers thyroid cancer susceptibility through the recruitment of USF1/USF2 transcription factors. PLoS Genet (2009) 5:e1000637. doi: 10.1371/journal.pgen.1000637

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Damiola F, Byrnes G, Moissonnier M, Pertesi M, Deltour I, Fillon A, et al. Contribution of ATM and FOXE1 (TTF2) to risk of papillary thyroid carcinoma in Belarusian children exposed to radiation. Int J Cancer (2014) 134:1659–68. doi: 10.1002/ijc.28483

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Penna-Martinez M, Epp F, Kahles H, Ramos-Lopez E, Hinsch N, Hansmann ML, et al. FOXE1 association with differentiated thyroid cancer and its progression. Thyroid (2014) 24:845–51. doi: 10.1089/thy.2013.0274

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Maillard S, Damiola F, Clero E, Pertesi M, Robinot N, Rachedi F, et al. Common variants at 9q22.33, 14q13.3, and ATM loci, and risk of differentiated thyroid cancer in the French Polynesian population. PLoS One (2015) 10:e0123700. doi: 10.1371/journal.pone.0123700

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zhuang Y, Wu W, Liu H, Shen W. Common genetic variants on FOXE1 contributes to thyroid cancer susceptibility: evidence based on 16 studies. Tumour Biol (2014) 35:6159–66. doi: 10.1007/s13277-014-1896-y

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Pereda CM, Lesueur F, Pertesi M, Robinot N, Lence-Anta JJ, Turcios S, et al. Common variants at the 9q22.33, 14q13.3 and ATM loci, and risk of differentiated thyroid cancer in the Cuban population. BMC Genet (2015) 16:22. doi: 10.1186/s12863-015-0180-5

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Matsuse M, Takahashi M, Mitsutake N, Nishihara E, Hirokawa M, Kawaguchi T, et al. The FOXE1 and NKX2-1 loci are associated with susceptibility to papillary thyroid carcinoma in the Japanese population. J Med Genet (2011) 48:645–8. doi: 10.1136/jmedgenet-2011-100063

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Tomaz RA, Sousa I, Silva JG, Santos C, Teixeira MR, Leite V, et al. FOXE1 polymorphisms are associated with familial and sporadic nonmedullary thyroid cancer susceptibility. Clin Endocrinol (Oxf) (2012) 77:926–33. doi: 10.1111/j.1365-2265.2012.04505.x

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Jones AM, Howarth KM, Martin L, Gorman M, Mihai R, Moss L, et al. Thyroid cancer susceptibility polymorphisms: confirmation of loci on chromosomes 9q22 and 14q13, validation of a recessive 8q24 locus and failure to replicate a locus on 5q24. J Med Genet (2012) 49:158–63. doi: 10.1136/jmedgenet-2011-100586

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Wang YL, Feng SH, Guo SC, Wei WJ, Li DS, Wang Y, et al. Confirmation of papillary thyroid cancer susceptibility loci identified by genome-wide association studies of chromosomes 14q13, 9q22, 2q35 and 8p12 in a Chinese population. J Med Genet (2013) 50:689–95. doi: 10.1136/jmedgenet-2013-101687

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Bonora E, Rizzato C, Diquigiovanni C, Oudot-Mellakh T, Campa D, Vargiolu M, et al. The FOXE1 locus is a major genetic determinant for familial nonmedullary thyroid carcinoma. Int J Cancer (2014) 134:2098–107. doi: 10.1002/ijc.28543

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Rogounovitch TI, Bychkov A, Takahashi M, Mitsutake N, Nakashima M, Nikitski AV, et al. The common genetic variant rs944289 on chromosome 14q13.3 associates with risk of both malignant and benign thyroid tumors in the Japanese population. Thyroid (2015) 25:333–40. doi: 10.1089/thy.2014.0431

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Liyanarachchi S, Wojcicka A, Li W, Czetwertynska M, Stachlewska E, Nagy R, et al. Cumulative risk impact of five genetic variants associated with papillary thyroid carcinoma. Thyroid (2013) 23:1532–40. doi: 10.1089/thy.2013.0102

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Saenko VA, Rogounovitch TI. Genetic Polymorphism Predisposing to Differentiated Thyroid Cancer: A Review of Major Findings of the Genome-Wide Association Studies. Endocrinol Metab (Seoul) (2018) 33:164–74. doi: 10.3803/EnM.2018.33.2.164

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Li J, Wang X, Dong J. Association of rs6983267 Polymorphism and Thyroid Cancer Susceptibility: A Systematic Review and Meta-Analysis. Med Sci Monit (2016) 22:1866–71. doi: 10.12659/msm.896507

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lloyd RV, Osamura RY, Kloppel G, Rosai J eds. WHO classification of tumours of endocrine organs. Lyon: IARC Press (2017).

Google Scholar

31. Sobin LH, Gospodarowicz MK, Wittekind C, International Union against Cancer. TNM classification of malignant tumours. 7th edition. Chichester: Wiley-Blackwell (2010).

Google Scholar

32. Chang CC, Chow CC, Tellier LCAM, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience (2015) 4:s13742-015-0047-8. doi: 10.1186/s13742-015-0047-8

CrossRef Full Text | Google Scholar

33. Jendrzejewski J, Liyanarachchi S, Nagy R, Senter L, Wakely PE, Thomas A, et al. Papillary Thyroid Carcinoma: Association Between Germline DNA Variant Markers and Clinical Parameters. Thyroid (2016) 26:1276–84. doi: 10.1089/thy.2015.0665

PubMed Abstract | CrossRef Full Text | Google Scholar

34. He H, Li W, Liyanarachchi S, Jendrzejewski J, Srinivas M, Davuluri RV, et al. Genetic predisposition to papillary thyroid carcinoma: involvement of FOXE1, TSHR, and a novel lincRNA gene, PTCSC2. J Clin Endocrinol Metab (2015) 100:E164–72. doi: 10.1210/jc.2014-2147

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Jendrzejewski J, He H, Radomska HS, Li W, Tomsic J, Liyanarachchi S, et al. The polymorphism rs944289 predisposes to papillary thyroid carcinoma through a large intergenic noncoding RNA gene of tumor suppressor type. Proc Natl Acad Sci U S A (2012) 109:8646–51. doi: 10.1073/pnas.1205654109

PubMed Abstract | CrossRef Full Text | Google Scholar

36. He H, Li W, Liyanarachchi S, Wang Y, Yu L, Genutis LK, et al. The role of NRG1 in the predisposition to papillary thyroid carcinoma. J Clin Endocrinol Metab (2018) 103:1369–79. doi: 10.1210/jc.2017-01798

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Graffelman J, Chang C, Puig X, Wigginton J. Statistical Tests and Graphics for Hardy-Weinberg Equilibrium (2019). Available at: https://cran.r-project.org/web/packages/HardyWeinberg/HardyWeinberg.pdf (Accessed March 12, 2020).

Google Scholar

38. Cassell DL. A Randomization-test Wrapper for SAS PROCs. Available at: https://support.sas.com/resources/papers/proceedings/proceedings/sugi27/p251-27.pdf (Accessed March 12, 2020).

Google Scholar

39. Nikitski A, Saenko V, Shimamura M, Nakashima M, Matsuse M, Suzuki K, et al. Targeted Foxe1 Overexpression in Mouse Thyroid Causes the Development of Multinodular Goiter But Does Not Promote Carcinogenesis. Endocrinology (2016) 157:2182–95. doi: 10.1210/en.2015-2066

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang F, Yan D, Ji X, Han J, Chen M, Qiao H, et al. rs965513 polymorphism as a common risk marker is associated with papillary thyroid cancer. Oncotarget (2016) 7:41336–45. doi: 10.18632/oncotarget.9324

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Nikitski AV, Rogounovitch TI, Bychkov A, Takahashi M, Yoshiura KI, Mitsutake N, et al. Genotype Analyses in the Japanese and Belarusian Populations Reveal Independent Effects of rs965513 and rs1867277 but Do Not Support the Role of FOXE1 Polyalanine Tract Length in Conferring Risk for Papillary Thyroid Carcinoma. Thyroid (2017) 27:224–35. doi: 10.1089/thy.2015.0541

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Montero-Conde C, Ruiz-Llorente S, Dominguez JM, Knauf JA, Viale A, Sherman EJ, et al. Relief of feedback inhibition of HER3 transcription by RAF and MEK inhibitors attenuates their antitumor effects in BRAF-mutant thyroid carcinomas. Cancer Discov (2013) 3:520–33. doi: 10.1158/2159-8290.CD-12-0531

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Bodmer D, Schepens M, Eleveld MJ, Schoenmakers EF, Geurts van Kessel A. Disruption of a novel gene, DIRC3, and expression of DIRC3-HSPBAP1 fusion transcripts in a case of familial renal cell cancer and t ()(q35;q21). Genes Chromosomes Cancer (2003) 38:107–16. doi: 10.1002/gcc.10243

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Coe EA, Tan JY, Shapiro M, Louphrasitthiphol P, Bassett AR, Marques AC, et al. The MITF-SOX10 regulated long non-coding RNA DIRC3 is a melanoma tumour suppressor. PLoS Genet (2019) 15:e1008501. doi: 10.1371/journal.pgen.1008501

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Michailidou K, Hall P, Gonzalez-Neira A, Ghoussaini M, Dennis J, Milne RL, et al. Large-scale genotyping identifies 41 new loci associated with breast cancer risk. Nat Genet (2013) 45:353–61, 361e1-2. doi: 10.1038/ng.2563

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Shen Z, Ren W, Bai Y, Chen Z, Li J, Li B, et al. DIRC3 and near NABP1 genetic polymorphisms are associated laryngeal squamous cell carcinoma patient survival. Oncotarget (2016) 7:79596–604. doi: 10.18632/oncotarget.12865

PubMed Abstract | CrossRef Full Text | Google Scholar

47. NCBI. Gene database. Available at: https://www.ncbi.nlm.nih.gov/gene/60625#gene-expression (Accessed March 12, 2020).

Google Scholar

48. Yi J, Zhou LY, Yi YY, Zhu X, Su XY, Zhao Q, et al. Low Expression of Pseudogene POU5F1B Affects Diagnosis and Prognosis in Acute Myeloid Leukemia (AML). Med Sci Monit (2019) 25:4952–9. doi: 10.12659/MSM.914352

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Pan Y, Zhan L, Chen L, Zhang H, Sun C, Xing C. POU5F1B promotes hepatocellular carcinoma proliferation by activating AKT. Biomed Pharmacother (2018) 100:374–80. doi: 10.1016/j.biopha.2018.02.023

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Shen L, Du M, Wang C, Gu D, Wang M, Zhang Q, et al. Clinical significance of POU5F1P1 rs10505477 polymorphism in Chinese gastric cancer patients receving cisplatin-based chemotherapy after surgical resection. Int J Mol Sci (2014) 15:12764–77. doi: 10.3390/ijms150712764

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Hayashi H, Arao T, Togashi Y, Kato H, Fujita Y, De Velasco MA, et al. The OCT4 pseudogene POU5F1B is amplified and promotes an aggressive phenotype in gastric cancer. Oncogene (2015) 34:199–208. doi: 10.1038/onc.2013.547

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Fu C, Xu X, Lu W, Nie L, Yin T, Wu D. Increased expression of long non-coding RNA CCAT2 predicts poorer prognosis in patients with hepatocellular carcinoma. Med (Baltimore) (2019) 98:e17412. doi: 10.1097/MD.0000000000017412

CrossRef Full Text | Google Scholar

53. Shah MY, Ferracin M, Pileczki V, Chen B, Redis R, Fabris L, et al. Cancer-associated rs6983267 SNP and its accompanying long noncoding RNA CCAT2 induce myeloid malignancies via unique SNP-specific RNA mutations. Genome Res (2018) 28:432–47. doi: 10.1101/gr.225128.117

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Yan L, Wu X, Yin X, Du F, Liu Y, Ding X. LncRNA CCAT2 promoted osteosarcoma cell proliferation and invasion. J Cell Mol Med (2018) 22:2592–9. doi: 10.1111/jcmm.13518

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Bai JG, Tang RF, Shang JF, Qi S, Yu GD, Sun C. Upregulation of long noncoding RNA CCAT2 indicates a poor prognosis and promotes proliferation and metastasis in intrahepatic cholangiocarcinoma. Mol Med Rep (2018) 17:5328–35. doi: 10.3892/mmr.2018.8518

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Van den Broeck T, Joniau S, Clinckemalie L, Helsen C, Prekovic S, Spans L, et al. The role of single nucleotide polymorphisms in predicting prostate cancer risk and therapeutic decision making. Biomed Res Int (2014) 2014:627510. doi: 10.1155/2014/627510

PubMed Abstract | CrossRef Full Text | Google Scholar