Multidimensional data analysis revealed thyroiditis-associated TCF19 SNP rs2073724 as a highly ranked protective variant in thyroid cancer

Genome-wide and phenome-wide associations between TCF19 and immune diseases and human cancers

We collected multiple GWAS studies, including UKB, FinnGen, and GWAS catalogue, and identified that TCF19 is significantly associated with multiple autoimmune diseases and human cancers, for example, type 1 diabetes and autoimmune thyroid diseases (P=1.0x10^-23), cervical cancer (4.0x10^-25) and lung cancer (5.0x10^-19) (Table 1). We also observed that rs2073724, a potential deleterious missense variation predicted by Sorting Intolerant From Tolerant (SIFT) [16], is associated with thyroid function (P=2.9x10-19), thyroiditis (P=2.3x10-7), hypothyroidism (P=1.0x10-8) and multiple autoimmune diseases (P=2.7x10^-10) (Table 2).

Moreover, the GTEx dataset showed that compared with the wild type allele (Pro), the rs2073724 alternative allele (Leu) increased the expression level of TCF19 in almost all human tissues (thyroid tissue, NES=0.99, P=9.6x10^-66, Supplementary Figure 1). The causal relationship between the alternative allele and increased aetiology of thyroiditis and thyroid cancer. Through analysis of the TCF19-related SNP in UCSC, we found that rs2073724 was located in exon 3 of TCF19 (Supplementary Figure 2). Hence, exonic SNPs are important SNPs, and the function of exonic SNPs in thyroid cancer remains largely unknown. Therefore, the potential role of TCF19 and its exonic SNP (rs2073724) in thyroid cancer tumorigenesis and the specific mechanism drew our attention.

TCF19 is upregulated in thyroid cancer and correlated with cancer progression and a poor prognosis

To further investigate the expression pattern of TCF19 in human cancers, we first analysed The Cancer Genome Atlas (TCGA) database and found that TCF19 mRNA levels were upregulated in most human cancers (meta-analysis, log2FC=0.68, P=2.84x10^-13) including thyroid cancer (Supplementary Figure 3 and Figure 1A).

Figure 1. High expression of TCF19 is associated with poor prognosis in thyroid cancer. (A) The expression pattern of TCF19 in total was analyzed in 512 THCA tissues and 59 normal controls (TCGA database). (B) Kaplan-Meier analysis evaluating the association between Progress Free Interval and TCF19 expression in thyroid cancer patients (N=512) in the TCGA-THCA cohort. Progress Free Interval analysis was further stratified by TCF19-High and TCF19-Low characteristics for Kaplan–Meier analysis. (C) The expression pattern of TCF19 in GSE29265 and GSE33630, comparing normal thyroid tissues with thyroid cancer. (D–G) The relationship between TCF19 expression and T stage (D), N stage (E), M stage (F) and TNM stage (G) in TCGA cohort. (H) The expression pattern of TCF19 in 20 THCA tumor and their adjacent tissues. *p<0.05, **p < 0.01, ***p<0.001, NS, not significant.

Furthermore, in two independent GEO datasets (GSE29265 and GSE33630), TCF19 mRNA levels were significantly higher in thyroid cancer than in matched normal controls (Figure 1B). To confirm the results of the datasets, we used a quantitative reverse transcription PCR (RT-PCR) method to quantify TCF19 mRNA levels in 10 pairs of thyroid cancer tissues and their corresponding adjacent normal thyroid tissues. Upregulation of TCF19 was detected in 10 pairs, accounting for 90% of all tested samples (Figure 1C). We also analyzed the TCF19 expression levels in relation to TNM stages in thyroid cancer patients using the TCGA database. We observed that TCF19 expression was elevated in T3 stage primary lesions as compared to T2 (Figure 1D). Similarly, patients with lymph node metastasis exhibited higher TCF19 levels than those without (Figure 1E). Although no significant differences were found in metastasis (M), this could be due to the small number of cases with distant metastasis (M1) in our sample (Figure 1F). Moreover, TCF19 expression was significantly higher in stage III and IV tumors compared to stage II (Figure 1G). Furthermore, patients with high TCF19 expression had shorter progression-free intervals (PFIs) than those with low TCF19 expression (HR=2.04, p=0.012) (Figure 1H). These results indicate that TCF19 is highly expressed in thyroid cancer and that patients with high TCF19 expression have a poor prognosis.

TCF19 promotes thyroid cancer progression, and the C>T variant of rs2073724 plays a protective role in thyroid cancer tumorigenesis

To detect the function of TCF19 and the SNP in thyroid cancer progression, we established TCF19[C] (rs2073724[C]) expression and TCF19[T] (rs2073724[T]) overexpression in thyroid cancer. RT-qPCR and Western blot analysis confirmed that TCF19 was strongly overexpressed. There was no significant difference in the expression level between TCF19[C] and TCF19[T] (Figure 2A and Supplementary Figure 4A). CCK-8 and EdU assays revealed that TCF19[C] overexpression promotes thyroid cancer cell proliferation, and TCF19[T] partly blocked these effects (Figure 2B, 2C). Additionally, we also performed Transwell and wound healing assays. As we observed previously TCF19[C] significantly increased the migration and invasion of thyroid cancer cells, while overexpression of TCF19[T] notably reduced these abilities affected by TCF19[C] overexpression (Figure 2D). Taken together, TCF19 plays an oncogenic role in thyroid cancer progression. Notably, the alternative allele of the SNP is required for the effect of TCF19 on thyroid cancer progression. In addition, we also established stable TCF19-knockdown TC cells (Supplementary Figure 4A, 4B). As expected, downregulation of TCF19 suppressed the proliferation, migration, and invasion of thyroid cancer cells compared with control cells, as demonstrated by CCK-8, EdU, Transwell, and wound healing assays (Supplementary Figure 4C–4F). These results indicate that TCF19 regulates the proliferation, migration and invasion of thyroid cancer cells in vitro. To investigate TCF19’s function in tumorigenesis in vivo, we conducted subcutaneous tumor formation assays using nude mice. We injected nude mice subcutaneously with C643 cells overexpressing either TCF19[C] or TCF19[T], along with a control group. The mice were euthanized 24 days post-injection to assess tumor development. Our findings revealed that TCF19[C] overexpression facilitated the growth of subcutaneous tumors. In contrast, TCF19[T] overexpression was observed to impede tumor progression when compared to TCF19[C] (Figure 2F–2H). More importantly, the C>T variant of rs2073724 diminishes the carcinogenic effect of TCF19.

Figure 2. TCF19[C] and TCF19[T] expression impacts thyroid cancer progression. (A) TCF19 expression was detected in KTC-1 and C643 cells transfected with pcDNA, TCF19[C] and TCF19[T] plasmids by Western blotting. (B) CCK8 assay was used to detect cell viability. (C) The EdU method was used to detect DNA synthesis. (D) Wound healing assay was used to detect cell migration. (E) Transwell assay was used to detect cell invasion and migration in the cells indicated above. (F) Tumor images in tumor-bearing mouse groups. (G) Mouse tumor weights in various groups. (H) Mouse tumor volumes were measured every four days. Data are shown as means ± S.D. *p < 0.05, **p < 0.01, ***p < 0.001, NS, not significant.

Identification of TCF19[C] and TCF19[T] targets by RNA-seq

To identify the potential targets regulated by TCF19, RNA sequencing (RNA-seq) was performed in control, TCF19[C] and TCF19[T] overexpression cells. As shown in Supplementary Figure 5A, PCA analysis indicates that the different RNA-seq samples could be completely distinguished. Overexpression of TCF19[C] resulted in abnormal gene expression, with 339 upregulated genes and 424 downregulated genes. However, when TCF19[T] was overexpressed, a total of 723 genes were upregulated and 1865 genes were downregulated in these cells compared to TCF19[C]-overexpressing cells (Figure 3A, 3B). The results indicate that TCF19[C] and TCF19[T] affect the expression of many genes in thyroid cancer cells. GO and KEGG analyses showed that the differentially expressed genes (DEGs) were mainly enriched in cancer-related programs or signalling pathways, such as the TNF signalling pathway, IL-17 signalling pathway and cytokine-cytokine receptor interaction (Figure 3C, 3D). And then we found that there were 31 genes associated with SNP (Figure 3E and Supplementary Figure 5B). Additionally, GO and KEGG analyses revealed enrichment of signatures associated with signalling receptor activator activity, the IL-17 signalling pathway, and the TNF signalling pathway (Figure 3F, 3G). Based on these results, TCF19 appears to be a potent regulator of thyroid cancer progression, especially in regard to the inflammatory response. In addition, the SNP rs2073724 plays a crucial role in maintaining protein function.

Figure 3. Identification of TCF19[C] and TCF19[T] targets by RNA-seq. (A) Volcano plots displaying differentially expressed genes (DEGs) in microarray data comparing PCDH with TCF19[C] and TCF19[C] with TCF19[T] KTC-1 cells. The numbers of significantly variant genes (FC > 2.0, P < 0.05) were shown. Vertical dashed lines indicate cut-off of FC (2.0), whereas the horizontal dashed lines indicate cut-off of p-value (0.05). (B) Heatmap of DEGs identified by RNA-seq. (C) Gene Ontology (GO) enrichment analysis of DEGs, p-value<0.05. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DEGs, p-value<0.05. (E) Heatmap analyzed from 31 total genes obtained by the intersection of cluster 1 (PCDH vs TCF19[C]) and cluster 2 (TCF19[C] VS TCF19[T]). (F) GO analysis of genes described in (E). (G) KEGG analysis of genes described in (E).

The tumour microenvironment (TME) is a complex and evolving system. Tumour occurrence, growth and metastasis strongly depend on the internal environment of tumour cells. It is related not only to the structure, function and metabolism of tumour tissues but also to the internal environment of tumour cells [17, 18]. To further stress the role of TCF19 in regulating the TME of thyroid cancer. We analysed the TCGA databases and found a positive correlation between TCF19 and immune scores (Supplementary Figure 6A) and stromal scores (Supplementary Figure 6B) and a negative correlation between TCF19 and tumour purity (Supplementary Figure 6C) in thyroid cancer. We further performed xCell and EPIC analyses to explore whether TCF19 expression correlates with immune cell infiltration. The results showed that THCA samples with high TCF19 expression had higher levels of monocyte and macrophage infiltration than samples with low TCF19 expression (Supplementary Figure 6D, 6E). Next, we quantified 8 immune cell populations and 2 nonimmune stromal populations (endothelial cells and fibroblasts) by MCP-counter scoring. The results showed that the proportions of B lineage, monocytic lineage, CD8+ T cells, cytotoxic lymphocytes, fibroblasts and myeloid dendritic cells were higher, while the proportion of endothelial cells was lower in the TCF19 high group than in the TCF19 low group (Supplementary Figure 6F). We also used the Tumour Immune Dysfunction and Exclusion (TIDE) algorithm [19] to predict the response to immune checkpoint blockade therapy (ICB). In general, high TIDE score was associated with a lower ICB response rate. Our results suggest that TCF19 overexpression may confer resistance to ICB therapy (Supplementary Figure 6G).

TCF19[T] reduces its DNA binding ability, thus affecting protein function

To further validate that the SNP mutation could abolish TCF19 protein function, we performed RT-PCR to detect the downstream genes of TCF19 from the RNA-seq. The results indicated that TCF19 could significantly regulate its target genes, while TCF19[T] partly abolished this function (Figure 4A). To further stress this interesting question, we analysed the SNP mutation site of the protein. We found that the SNP mutation site was located in the proline-rich region of TCF19 (Supplementary Figure 7A, 7B). A previous study indicated that TCF19 contains a proline-rich region, which is a common characteristic of transactivating factors [20]. We noticed that the SNP was located in the proline-rich region of the TCF19 protein. Therefore, we hypothesized that the SNP mutation affects the transactivating ability of TCF19. Then, we analysed our unpublished CUT&Tag data for TCF19 and found that TCF19 could directly bind to the promoter of the SRSF6 gene (Figure 4B). To further verify whether TCF19 activates the SRSF6 promoter through direct binding, we performed chromatin immunoprecipitation (ChIP) assays. As expected, TCF19 was recruited to the promoter sites. However, TCF19[T] strongly reduced the binding ability at the promoter of SRSF6 (Figure 4C). Taken together, TCF19 may be an exciting new therapeutic target for aggressive thyroid cancer. The alternative allele of the rs2073724 variant disrupts the oncogenic functions of TCF19, implicating the rs2073724 C>T variant as a causal genetic variant influencing thyroid cancer risk.

Figure 4. TCF19[T] reduces its DNA binding ability, thus affecting protein function. (A) CSF3, SRSF6, NCOA5, IL1B, NFE2, MAT2A, TXNIP, GALNT3 mRNA expression levels in cells stably transfected with PCDH, TCF19[C] and TCF19[T] plasmids were measured by qRT-PCR, respectively. (B) CUT&Tag data showed TCF19 could directly bind to the promoter of SRSF6 gene. (C) ChIP-qPCR analysis revealed potential TCF19-binding sites in the SRSF6 promoter region. After SNP mutation, the binding ability was weakened.

Bioinformatics analysis

Previous GWAS and pheWAS studies were collected from multiple sources, including the GWAS catalogue, Opentargets and ExPheWas. The GTEx’s final dataset (V8) contains DNA data from 838 postmortem donors and 17,382 RNA-seq across 54 tissue sites and two cell lines where eQTLs have been well identified and used in previous research. Cancer samples were collected from the TCGA project (N = 10,490). The gene expression level was log2 transformed before the meta-analysis. A fixed effect model and a random effect model were both applied for the aggregation. The 95% CI was applied to show the risk and protective effect on overall survival time. To show more details for different studies, any standardized mean difference (SMD) higher than 3 and lower than -3 is shown with an arrow. Blue filled parallelograms represent the SMD for the fixed effect model and random effect model.

Cell lines and cell culture

KTC-1 and C643 cells were identified by STR analysis. All cell lines were cultured in RPMI-1640 or DMEM (Gibco, USA) supplemented with 10% foetal bovine serum, penicillin (Gibco, USA) and streptomycin (Gibco) and kept in a humid atmosphere of 5% CO₂ at 37° C.

Plasmid construction and lentivirus infection

The cDNA of the TCF19 CDS region was amplified and cloned and inserted into the lentiviral vector pCD-Puro-3×Flag. Catalytic mutant TCF19 was synthesized by GENEWIZ using the same vector as wild-type TCF19. The shRNA targeting human TCF19 was cloned and inserted into pLV-shRNA-puro-mCherry. All constructed vectors were verified by DNA sequencing. The target sequences were as follows: shTCF19-1, 5’-CCAAGGUACUUUGGUCAAUAA-3’; shTCF19-2, 5’-CCCAGGAAGAAACUCCGUGUA-3’; and shTCF19-3, 5’-ACACUGAUCCUAAACUCCAUA-3’. Plasmids were transfected into 293T cells using PEI, Opti-MEM I reduced serum medium (Gibco, USA) and packaging vectors (PAX8 and PVSVG). Forty-eight hours and 72 hours after transfection, viral supernatants were collected, passed through a 0.45 μm filter, and infected at 50% confluency to target cells.

RNA extraction and RT-qPCR

Total RNA was isolated using RNA Extraction Solution (Servicebio, Wuhan, China), followed by cDNA synthesis using HiScript III RT SuperMix for qPCR (+g DNA wiper) (Vazyme Biotech, Nanjing, China). RNA expression levels were determined by SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China) on a Bio-Rad QX100 Droplet Digital PCR System (USA), and the RNA expression of each target was calculated by the 2^-ΔΔCT method and normalized to that of β-actin. qPCR was performed as shown in Supplementary Table 1.

Western blotting

Total protein was extracted with RIPA lysis buffer (Solarbio, Beijing, China) containing protease inhibitors. Proteins were detected using BCA protein detection kit (Biosharp, Beijing, China), then separated by SDS-PAGE and migrated to PVDF membrane (Millipore, USA). The membrane was then incubated with primary antibody against GAPDH (1:5000) (Boster, USA, A00227-1), TCF19 (1:1000) (Abcam, UK, ab230005), and FLAG (1:1000) (Sigma, USA, F1804). Western blot bands were detected using a Sparkjade ECL super (Sparkjade, Shandong, China) and an imaging system (Bio-Rad, USA).

CCK-8 assay

THCA cells were seeded in 96-well plates at a density of 1 × 10³ cells per well and then treated with CCK-8 reagent (APExBIO, USA). The absorbance was measured at 450 nm.

EdU assay

Cells were cultured in 24-well plates and treated with 1 ml medium containing 1 μl EdU (10 mM). After incubation for 2 hours at 37° C and 5% CO₂, the cells were fixed with 4% paraformaldehyde for 30 minutes and incubated with 0.5% Triton-X-100 in PBS for 20 minutes. The nuclei were stained with Hoechst 33342 dye. The proliferation rate was calculated according to the manufacturer’s instructions (Bioscience, Shanghai, China). In each group, three regions were randomly selected and imaged with a fluorescence microscope (Leica, Wetzlar, Germany).

Wound healing assay

Cells were seeded and cultured until a 90% confluent monolayer formed and were then scraped with sterile pipette tips and processed in FBS-free medium as indicated in the text. Three areas were randomly selected under the microscope at 0 h and 48 h, and the distance of the cells that migrated to the scratched area was measured.

Migration and invasion assays

Migration or invasion assays were performed using 24-well plates immersed in 8.0 μm pore size Transwell filters (Corning, USA) with or without precoated diluted Matrigel (Corning, USA). THCA cells (2×10^4) in serum-free medium were added to the upper chamber, and medium containing 10% foetal bovine serum was added to the lower chamber. After incubation at 37° C for 24 h (migration) or 48 h (invasion), cells were fixed on slides and stained with crystal violet (Solarbio, Beijing, China). Infiltrating cells in 3 random fields were then counted under the microscope.

RNA-seq data analysis

Total RNA was harvested and isolated from stable TCF19 overexpression, TCF19 mutation and scramble control cells. The library was sequenced on an Illumina NovaSeq 6000. Differentially expressed genes were defined as those with fold change >2 or <0.5, p<0.05 and were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

ChIP-PCR assay

ChIP detection was performed using the SimpleChIP®Plus Enzymatic Chromatin IP Kit (magnetic beads) (CST, USA). For ChIP-qPCR detection, the primers designed for the promoter region of SRSF6 are listed in Supplementary Table 1.

Statistical analysis

All the data presented are from at least three independent biological experiments and are shown as the mean ± SD. An unpaired two-tailed Student’s t test was used for statistical analysis of the experiments. P<0.05 was considered to indicate statistical significance. All data were statistically analysed with GraphPad Prism 7.0.