Papillary Thyroid Cancer Remodels the Genomic Fabrics of DNA Replication, Repair and Transcription

1. Introduction

The American Cancer Society estimates that in 2024 U.S.A. will register new 44,020 (31,520 women and 12,500 men) cases of thyroid cancer (TC), out of which, and 2,170 (1,180 females + 990 males) might die because of it [1]. The most lethal form is the non-differentiated (anaplastic) TC with 5 months’ median survival rate [2]. The least lethal is the differentiated papillary thyroid carcinoma (PTC) that covers over 8 out of 10 thyroid malignancies, growths slowly and has a favorable (> 95%) 5-year survival prognosis [3,4]. However, not treated in time, PTC can invade adjacent structures and metastasize in neck lymph nodes [5,6].

It is considered that the T1799A point mutation of the BRAF gene is responsible for 45% of the PTC cases [7] albeit mutations of several other genes (e.g.: RET, NTRK1, TP53) were also reported e.g.: [8,9,10,11] as present in numerous PTC cases. In total, the 39.0 Release of the National Cancer Institute GDC Data Portal (Dec. 04, 2023) lists 5,869 mutations affecting 19,313 genes in the registered 1,121 cases of thyroid cancer [12]. In addition to specific DNA mutations, regulation of the expression levels of selected gene sets were also pointed out as PTC transcriptomic signatures (e.g. [13,14,15]).

In order to quantify the power of the expression a particular gene to regulate the cancer phenotype, we have introduced the Gene Commanding Height (GCH) score [16]. Nonetheless, our PTC transcriptomic study [17] revealed that 1/3 of 544 known cancer biomarkers had higher GCH in the malignant region of the tumor, 1/3 in the benign region, while the rest did not discriminate between the two regions. This result raised the general question of the biomarkers’ utility for cancer diagnostic [18] and therapy [19], especially because, together with the blamed biomarker(s), sequence and/or expression of numerous other protein coding (e.g. [20]) and non-coding (e.g. [21]) genes are altered in cancer [22]. Moreover, all our cancer genomics studies revealed that biomarkers have low GCH scores and, as minor players in the cell life, their manipulation might be of little consequences for the survival and proliferation of the malignancy.

Here, we report the PTC-induced changes in the KEGG-constructed functional pathways [23] responsible for the genetic information processing in a surgically removed PTC tumor. The study complements the not so rich literature related to the genomic alterations in the DNA replication [24], repair [25,26,27,28] and transcription [29,30], chromatin remodeling [31], translation, and protein processing [32] and trafficking [33].

Gene expressions in the cancer nodule and surrounding normal thyroid tissue were compared with the transcriptomic profiles of the papillary (BCPAP, [34]) and anaplastic (8505C, [35]) human thyroid cancer cell cultures. As already adopted by us [36] and several other laboratories (e.g. [37,38,39]), the surrounding normal tissue in the tumor is the best reference to understand the genomic changes in the cancer nodule(s).

The analyses were carried out from the Genomic Fabric Paradigm (GFP) perspective that provides the most theoretically possible characterization of the transcriptome [40] and is able to identify the most legitimate targets for personalized anti-cancer gene therapy [41]. In addition to the traditional average expression level (AVE), GFP characterizes each gene by the Relative Expression Variation (REV) across biological replicas and Expression Coordination (COR) with each other gene in the same condition.

By comparing AVEs in the cancer nodule and the surrounding normal thyroid tissue one determines what gene was significantly up-/down-regulated according to certain cut-off criteria for the absolute fold-change and p-value of the t-test of means equality. AVE analysis is used also to identify the turned on/off genes.

REV analysis identifies the genes whose expression levels are the most/least controlled by the homeostatic mechanisms. Indirectly, REV points out the cell priorities in keeping within narrow intervals the fluctuations of the expression levels of critical genes to preserve the cell phenotype caused by nonhomogeneous and variable local environmental factors.

Moreover, COR analysis based on the Principle of Transcriptomic Stoichiometry [42]

stating that genes whose encoded products interact are coordinately expressed to

maximize the pathway efficiency. Thus, COR analysis singles out the most probable, statistically (p < 0.95) significant gene networking in functional pathways.

2. Materials and Methods

2.1. Gene Expression Data

Transcriptomic raw data obtained through our previous microarray experiments on thyroid cancer samples were downloaded from the publicly available Gene Expression Omnibus (GEO) of the National Center for Biotechnology Information (NCBI). As described in [17], 4 small pieces of malignant (here after denoted as T1, T2, T3, T4) and 4 of non-malignant (denoted as N1, N2, N3, N4) regions were collected from a surgically removed 32.0 mm pathological stage pT3NOMx [43] PTC tumor of a 33y old Asian woman. Tumor gene expression data [44] were compared with those from 4 papillary BCPAP (denoted as: Φ1, Φ2, Φ3, Φ4) and 4 anaplastic 8505C (denoted as: Θ1, Θ2, Θ3, Θ4) thyroid cancer cell cultures [45]. All thyroid cancer samples were profiled by us using Agilent Agilent-026652 Whole Human Genome Microarray 4x44K v2 [46]. The wet protocol and the raw data are fully described in the above cited, publicly accessible, GEO deposits.

2.2. Transcriptomic Analyses

The microarray data were filtered (eliminated all spots with foreground fluorescence less than twice the background in one microarray), normalized (to the median of valid spots) and analyzed using our standard GFP algorithms (presented in [47]). Besides the primary characteristics of adequately profiled individual genes in each condition: AVE, REV and COR (Appendix A), we considered also the derived characteristics: Relative Expression Control (REC) and Coordination Degree (COORD) (Appendix B).

The AVE, REV, COR, REC and COORD alterations in the malignant region of the tumor were quantified according to the algorithms from Appendix C and averaged for the genes included in the selected KEGG-constructed functional pathways. The arbitrarily introduced absolute fold-change cut-off (e.g. 1.5x) to consider a gene as significantly regulated might be too strict for very stably expressed genes or two laxed for the very unstably expressed ones across biological replicas. Therefore, we computed the absolute fold-change cut-off for each quantified gene considering the combined contributions of the biological variability in the compared conditions and the technical noise of the probing microarray spots [48].

Nevertheless, the transcriptomic distance (i.e. the Euclidian distance in the 3D orthogonal space of transcriptomic changes (WIR, ΔREC, ΔCOORD) with respect to the normal tissue (defined in C6) is the most comprehensive measure of the cancer-induce transcriptomic alteration of the thyroid.

2.3. Functional Pathways

We analyzed the following KEGG-constructed functional pathways responsible for the Genetic Information Processing in homo sapiens (hsa):

2.3.1. Transcription

(POL) hsa03020 RNA polymerase [49]

(BTF) hsa03022 Basal transcription factors [50]

(SPL) hsa03040 Spliceosome [51]

2.3.2. Translation

(RIB) hsa03010 Ribosome [52]

(AMI) hsa00970 Aminoacyl-tRNA biosynthesis [53]

(NCT) hsa03013 Nucleocytoplasmic transport [54]

(SUR) hsa03015 mRNA surveillance pathway [55]

(RBE) hsa03008 Ribosome biogenesis in eukaryotes [56]

2.3.3. Folding, Sorting and Degradation

(PEX) hsa03060 Protein export [57]

(PPE) hsa04141 Protein processing in endoplasmic reticulum [58]

(SIV) hsa04130 SNARE interactions in vesicular transport [59]

(UMP) hsa04120 Ubiquitin mediated proteolysis [60]

(SRS) hsa04122 Sulfur relay system [61]

(PRO) hsa03050 Proteasome [62]

(RND) hsa03018 RNA degradation [63]

2.3.4. Replication and Repair

(DER) hsa03030 DNA replication [64]

(BER) hsa03410 Base excision repair [65]

(NER) hsa03420 Nucleotide excision repair [66]

(MIR) hsa03430 Mismatch repair [67]

(HOR) hsa03440 Homologous recombination [68]

(NHJ) hsa03450 Non-homologous end-joining [69]

(FAP) 03460 Fanconi anemia pathway [70]

2.3.5. Chromosome

(ACM) hsa03082 ATP-dependent chromatin remodeling [71]

(PRC) hsa03083 Polycomb repressive complex [72].

In addition, we have analyzed in the surgically removed PTC tumor the regulation of the genes modules responsible for the cancer cells’ survival and proliferation included in KEGG-constructed hsa05200Pathways in cancer (PAC, [73]).

3. Results

3.1. Independence of the Three Transcriptomic Characteristics of Individual Genes

Figure 1 illustrates the independence of the primary characteristics AVE, REV and COR of the 30 quantified DER genes in the normal thyroid tissue (N), papillary thyroid cancer tissue (T), papillary thyroid cancer cell line BCPAP (Φ) and anaplastic thyroid cancer cell line 8505C (Θ). Correlations with PCNA (proliferating cell nuclear antigen) were selected for the PCNA role in ensuring the rate and accuracy of DNA replication, damage repair, chromatin formation, and segregation of the sister chromatids [74]. The characteristics were determined with the definition algorithms described in Appendix A.

Using the three characteristics (whose independence is visually evident) increases by almost four orders of magnitude the amount of workable transcriptomic information that can be derived from any microarray (or RNA-sequencing) experiment. Thus, by quantifying 14,904 unigenes in this experiment, we got for each type of sample (N, T, Φ, Θ) 14,904 AVEs + 14,904 REVs + 111,057,156 CORs = 111,086,964 (i.e. 7,453.5x more values to be analyzed than the 14,904 AVEs).

It is interesting to note that, with respect to the normal thyroid tissue N, the cancer increased the average expression of these genes by 46% in T, 142% in Φ, and 146% in Θ. RNASEH2A (ribonuclease H2, subunit A), known for its role in protecting the genomic integrity and progression of prostate cancer [75], had the most remarkable increase in the cultured BCPAP (39x) and 8505C (49x) cancer cells.

On the other hand, the large, 91% REV of RNASEH2B (ribonuclease H2, subunit B) in the anaplastic cells shows the high expression addaptability of this gene to ensure the proper DNA replication regardless of the tumor heterogeneous cell subpopulations.

Correlation analysis revealed that within the DER pathway, PCNA has 5 statistically significant synergistically expressed partners in N, but 11 in T, 21 in Φ, and 16 in Θ. These results indicate substantial increase of synchronious expressions DNA replication genes in thyroid cancer. The numbers of antagonistically expressed partners are: 0 in N, 2 in T, 3 in Φ, and 0 in Θ.

3.2. Cancer-Induced Regulation of Gene Expression Profile

3.2.1. Measures of Regulation of Expression Level

Figure 2 presents three ways to quantify the cancer-related regulation of the expression level of 51 randomly selected out of 74 quantified genes encoding polymerases based on the algorithms described in Appendix C.

The traditional measure of the percentages of significantly up- and down-regulated genes (according to a more or less arbitrarily introduced criterion) considers only the regulated genes and as uniform +1 or -1 contributors to the transcriptomic alteration. Instead, the other two measures from Figure 2 apply to all genes while discriminating also their contributions to the overall transcriptomic regulation. Thus, POLR1C (polymerase (RNA) I polypeptide C, 30kDa) had the largest up-regulation in both cell cultures lines even its expression level was not significantly changed in the malignant part of the surgically removed tumor (x_T = -1.09; x_Φ = 198.71; x_Θ = 221.73). POLR1C was reported as one of the ten most up-regulated proteins when infiltrative gastric cancer regions were compared to the adjacent normal tissue [76]. TAF1A (TATA-box binding protein associated factor, RNA polymerase I subunit A) was also a highly up-regulated gene in both cell lines although it was significantly down-regulated in the tumor (x_T = -2.12; x_Φ = 105.70; x_Θ = 128.88).

From the WIR perspective (Figure 2(c)), the largest negative contribution to the transcriptomic alteration of the polymerases were given in the two cell lines by PTRF (polymerase I and transcript release factor), although it appeared as not being significantly regulated in the tumor. The difference between high fold-change but low WIR for POLR1C and low absolute fold-change but high absolute WIR for PTRF in both cell lines comes from their difference in the average expression levels in the normal tissue: AVE^(POLR1C) = 0.17, AVE^(PTRF) = 16.80. (AVE is measured in the expression level of the median gene for that condition. POLR2L (RNA polymerase II, I and III subunit L) is another polymerase gene with very large negative contributions to the BCPAP and 8505C transcriptomes owing to the large expression level (AVE = 18.11) in the normal tissue.

3.2.2. Regulation of the KEGG-Constructed Functional Pathways Responsible for the Genetic Information Processing in the Malignat Region of the Thyroid Tumor

The regulated genes were identified using the cut-off criteria defined in Appendix C (C1) and the software “#PATHWAY#” described in Ref. [41].

Table 1 lists all the significantly up-regulated and down-regulated genes from the KEGG-constructed functional pathways responsible for Transcription and Translation in the cancer nodule of the surgically removed thyroid tumor.

Table 2 lists all the significantly up-regulated (bold) and down-regulated genes from the KEGG-constructed functional pathways responsible for Folding, Sorting and Degradation, Replication and Repair, and Chromosome in the cancer nodule of the surgically removed thyroid tumor.

3.2.3. Regulation of the Protein Processing in Endoplasmic Reticulum pathway

Figure 3 presents the localization of the regulated genes from the hsa04141 KEGG-constructed PPE pathway Protein Processing in Endoplasmic Reticulum [58] in the surgically removed PTC tumor. Within this gene subset, CRYAB (crystallin, alpha B), whose silencing might occur at the end of a stepwise dedifferentiation process in the thyroid gland [77], had the largest negative contribution in both PTC types of samples (x_T = -7.16, WIR_T = -105; x_Φ = -36.62, WIR_Φ = -610).

3.2.4. Regulation of the Cancer Cells’ Survival and Proliferation Genes

Figure 4 ^{modified frm} [73]: Localization of the regulated genes responsible for the cancer cells survival and proliferation within the KEGG-constructed functional pathways in cancer for the surgically removed PTC tumor. Regulated genes: BAX, BCL2, BIRC5 (baculoviral IAP repeat containing 5), CCND1 (cyclin D1), CDK6 (cyclin-dependent kinase 6), CDKN2A (cyclin-dependent kinase inhibitor 1A), CEBPA (CCAAT/enhancer binding protein (C/EBP), alpha), CKS2 (CDC28 protein kinase regulatory subunit 2), CSF1R (colony stimulating factor 1 receptor), DDB2 (damage-specific DNA binding protein 2), ESR1 (estrogen receptor 1), ETS1 (v-crk avian sarcoma virus CT10 oncogene homolog-like), FOS (FBJ murine osteosarcoma viral oncogene homolog), GTSK1/M1/M3/M5/o1/O2P1/T2 (glutathione S-transferase kappa 1/mu 1/mu 3/mu 5/omega 1/omega 2/pi 1/theta 2), HDAC1 (histone deacetylase 1), HEYL (hes-related family bHLH transcription factor with YRPW motif-like), HMOX1 (heme oxygenase (decycling) 1), JUN (jun proto-oncogene), MMP2 (matrix metallopeptidase 2 (gelatinase A, 72kDa gelatinase, 72kDa type IV collagenase)), MSH6 (mutS homolog 3), PAX8 (paired box 8), PDGF (platelet-derived growth factor beta polypeptide), PGF (placental growth factor), PIM1 (pim-1 oncogene), PPARG (peroxisome proliferator-activated receptor gamma), RUNX1T1 (runt-related transcription factor 1; translocated to, 1 (cyclin D-related)), RPS6KA5 (ribosomal protein S6 kinase, 90kDa, polypeptide 5), SLC2A1 (solute carrier family 2 (facilitated glucose transporter), member 1), SMAD3 (SMAD family member 3), SP1 (Sp1 transcription facto), STAT1/2/6 (signal transducer and activator of transcription 1), TCF7 (transcription factor 7 (T-cell specific, HMG-box)), TCF7L1 (transcription factor 7-like 1 (T-cell specific, HMG-box)), TGFB1 (transforming growth factor, beta 1), WNT11/3/3A/4/5A (wingless-type MMTV integration site family, member 11/3/3A/4/5A), ZBTB16 (zinc finger and BTB domain containing 17).

Genes in Figure 4 were identified by KEGG Pathways in Cancer [73] as responsible for blocking the differentiation, evading apoptosis, ensuring immortality and insensitivity to anti-growth signals, resistance to chemotherapy, as well proliferation and sustained angiogenesis. There is some degree of similarity between results in Figure 4 and those comparing gene expression profiles in the primary cancer nodule and the surrounding normal tissue from a surgically removed prostate tumor [78, Figure S1(a)]. For instance, GSTO2 and PGF were down-regulated while CSF1R, DDB2, HMOX1 and PGFB were up-regulated in both cancers. However, GSTP1 and PPARG were down-regulated in PTC but up-regulated in the profiled prostate tumor.

3.3. Additional Measures of Transcriptomic Alterations

Figure 5 illustrates our original novel measures of the cancer-induced transcriptomic alterations in expression control (ΔREC) and coordination degree (ΔCOORD), as well as the most comprehensive measure, the Transcriptomic Distance (TD). The measures, determined using the algorithms from Appendix C, were applied to the quantified 51 RNA polymerase II genes and their binding partners in the normal thyroid tissue (N), papillary thyroid cancer tissue (T), papillary thyroid cancer cell line BCPAP (Φ) and anaplastic thyroid cancer cell line 8505C (Θ).

Within this gene subset, the expression control (Figure 5(a)) had the largest increase for TAF7 (TATA-box binding protein associated factor 7) in T (ΔREC = 126), TAF6 (TATA-box binding protein associated factor 6) in Φ (ΔREC = 218) and POLR2B (RNA polymerase II subunit B) in Θ (ΔREC = 485). By contrast, the expression controls of POLR2L in T (ΔREC = -144) and Θ (ΔREC = -198), and GTF2E1 (general transcription factor IIE subunit 1) in Φ (ΔREC = -166) were substantially diminished. Interestingly, while the expression control of TAF6 increased in Φ, it stayed practically the same in T (ΔREC = 5), but decreased in Θ (ΔREC = -46), indicating significant changes in cancer cells’ priorities for controlling the expression fluctuations of this gene. The increased expression control of TAF7 is justified by its essential role for transcription, proliferation and differentiation [79].

The coordination degree (Figure 5(b)) was also largely affected by cancer. The largest increase was exhibited by ERCC3 (excision repair 3 TFIIH core complex helicase subunit) in T (ΔCOORD = 34) and Θ (ΔCOORD = 64), and TAF1L (TATA-box binding protein associated factor 1 like) in Φ (ΔCOORD = 60). Expression coordination of ERCC3, a.k.a. xeroderma-pigmentosum B (XPB) [80], increased also in Φ (ΔCOORD = 32). In contrast, GTF2I (general transcription factor IIi) in T (ΔCOORD = - 44), GTF2H4 (general transcription factor IIH subunit 4) in Φ (ΔCOORD = - 46) and TAF6 (TATA-box binding protein associated factor 6) in Θ (ΔCOORD = - 44) exhibited the largest decrease.

By the most comprehensive measure (Transcriptomic Distance, TD, (Figure 5(c)) POLR2L had the largest contribution to the overall cancer-induced transcriptomic alteration within this gene subset in T (TD^(T) = 145, TD⁽Φ⁾ = 497, TD^(Θ) = 495) and GTF2I in both cell lines (TD^(T) = 45, TD⁽Φ⁾ = 964 TD^(Θ) = 872). The least affected genes were: POLR2H (RNA polymerase II subunit H) in T (TD^(T) = 3, TD⁽Φ⁾ = 106, TD^(Θ) = 157), TBP (TATA-box binding protein) in Φ (TD^(T) = 15, TD⁽Φ⁾ = 26, TD^(Θ) = 18) and POLR2D (general transcription factor IIH subunit 1) in Θ (TD^(T) = 15, TD⁽Φ⁾ = 106, TD^(Θ) = 157). The average TD of this gene subset was 40 for the T samples, 99 for the BCPAP cells and 94 for the 8505C cells. It is important to note the substantially different contribution levels to the overall transcriptomic change of this gene subset in the three types of samples.

3.4. Cancer-Induced Remodeling of DNA Replication (DER) Pathway

Figure 6 presents the statistically (p < 0.95) significant synergism (COR > 0.950), antagonsim (COR < -0.950) and independence (|COR| < 0.050) of the expressions of genes involved in the sequentially connected multistep leading and lagging strands of the DNA replication [63].

3.5. Remodeling of the Coupling of Polymerase II Genes with Basal Transcription Factors

Figure 7 presents the statistically (p < 0.95) significant synergism, antagonsim and independence of the expressions of the polymerase II genes with their binding partners.

4. Discussion

In this report, we used the GFP [40] algorithms to reveal novel characteristics of the transcriptome that are ignored by the traditional analysis limited to the genes’ expression levels (e.g. [81,82]). Adding the independent Relative Expression Variation (REV) and Expression Correlation (COR) to the Average Expression Level (AVE) of individual genes increased by almost 75 hundred times the information obtained from our gene expression experiments on thyroid cancer samples [44,45]. The independence of the three characteristics, illustrated here for 30 DNA replication genes in all four types of samples, was confirmed for several other gene sets profiled on surgically removed tumors from kidney [47,48], prostate [36,78] and thyroid [22]. AVEs, REVs and CORs of all genes characterize the transcriptome as well as the numbers (AVE) of electronic devices of each type, their wiring (COR) and applied voltage fluctuations (REV) characterize a supercomputer.

REV (Appendix A3) and COR (Appendix A4) were used to define the derived characteristics of individual genes termed Relative Expression Control (REC, Appendix B1) and Coordination Degree (COORD, Appendix B2). The derived characteristics were used to define the novel measures of transcriptomic alteration termed Regulation of Expression Control (ΔREC; Appendix C4) and Regulation of the Coordination Degree (ΔCOORD; Appendix C5).

It is important to specify that the arbitrarily introduced (1.5x or 2.0x) absolute fold-change cut-off to consider a gene as significantly regulated was replaced by the flexible CUT (Appendix C1). CUT is computed for each transcript by evaluating the biological variability of the expression level across biological replicas and the technical noise of the probing spots in the compared conditions. Also important is that the expression levels of adequately quantified genes were normalized to the median expression of all genes in that condition which improves the accuracy of the conditions’ comparison.

Instead of the traditional percentages of up- and down-regulated genes that is limited to the significantly regulated, implicitly considering every affected gene as a uniform +1/-1 contributor to the expression profile alteration, we use the Weighted Individual (gene) Regulation (WIR; Appendix C2). WIR takes into account the total change of the expression level and the statistical confidence of the regulation. Altogether WIR, ΔREC and ΔCOORD were incorporated into the Transcriptomic Distance (TD; Appendix C6), the most theoretically possible comprehensive measure of the cancer-related transcriptomic change.

In the presented analyses, the transcriptomic regulations in the cancer nodule (T) and standard thyroid cancer cell lines BPAP (Φ) and 8505C (Θ) were referred to the gene expression profile of the surrounding normal tissue (N) in the surgically removed PTC tumor. While such comparison is natural for the cancer nodule, it might be disputable for the cell lines because of the differences in the cellular environment. In separate studies we proved the profound remodeling of the transcriptomes of each of two cell types (mouse cortical astrocytes and immortalized precursor oligodendrocytes) when co-cultured in insert systems [83,84]. These results suggest that a homo-cellular culture does not exactly repeat the genomic properties of the main cell type in the tissue.

The functional pathways analyzed in Table 1 and Table 2 are not mutually exclusive. For instance, POLD4 (polymerase (DNA-directed), delta 4, accessory subunit) is included in: hsa03030 DNA replication, hsa03410 Base excision repair, hsa03420 Nucleotide excision repair, hsa03430 Mismatch repair, and hsa03440 Homologous recombination. Interestingly, while we found POLD4 as significantly up-regulated (x = 1.90) in the malignant part of the tumor compared to the normal region, it was massively down regulated in both BPAP (x = -26.96) and X505C (x = -39.64) cells. A recent meta-analysis of TCGA (The Cancer Genome Atlas) data [85] has shown that POLD4 was significantly overexpressed in 17 types of cancer compared to adjacent normal tissue. Our results indicate again the limitations of the standardized cancer cell lines to replicate the characteristics of the real tumors owing to potential in vitro evolution and genomic instability caused by repeated passages of the cultured cells [86].

Figure 3 shows how the expressions of genes involved in the protein processing in endoplasmic reticulum were regulated in the cancerous part of the thyroid tumor. Interestingly, the ubiquitination gene UBXN1 (UBX domain protein 1), an important negative regulator of the unfolded protein response [87], was significantly up-regulated in the profiled tumor (x_T = 1.75) but massively down-regulated in both cell lines (x^(Φ) = -41.96, x^(Θ) = -106.57). While depletion of UBXN1 might be profitable for cancer cell survival and proliferation through protecting against endoplasmic reticulum (ER) stress [88,89] its up-regulation might facilitate migration and invasion of the cancer cells within the thyroid tissue as reported in prostate cancer [90]. Therefore, we speculate that, although in apparent contradiction, both UBXN1 up-regulation in the tissue (Figure 3) and UBXN1 down-regulation in the homo-cellular BPAP and 8505C cultures are justified by the molecular mechanisms controlled by the encoded protein. PTRF and UBE2O (ubiquitin-conjugating enzyme E2O) are other examples of genes that were not affected in the tumor but significantly regulated in both BPAP (x_PTRF = -31.50, x_UBE2O = 28.86) and 8505C (x_PTRF = -58.42, x_UBE2O = 11.70) cell lines. UBE2O was reported to ubiquitinate PTRF and inhibit the “secretion of exosome-related PTRF/CAVIN1” [91].

Figure 4 justifies the increased proliferation of the cancer cells through up-regulated CCND1 (x^(T) = 3.48; x^(Φ) = 4.78; x^(Θ) = 8.85), one of the 20 hub genes used as biomarker for the thyroid cancer [92]. Survival of the cancer cells is optimized through the up-regulation of several key genes ensuring evading apoptosis such are BIRC5 (x^(T) = 6.04) and DDB2 (x^(T) = 2.37). BIRC5 (x^(Φ) = 10.70; x^(Θ) = 7.35) is known for promoting multidrug resistance to chemotherapy [93] and the oncogenic roles of DDB2 were reported in a recent review [94]. However, with respect to the normal thyroid tissue, DDB2 was significantly downregulated in both cell lines (x^(Φ) = -3.01; x^(Θ) = -22.00). How to explain these opposite results? Within a hetero-cellular tissue (composed of normal and cancer thyrocytes and potential presence of immune cells), over-expression of genes such is DDB2 accelerates cancer cells proliferation and invasion of the thyroid. Nevertheless, expression of genes that are no longer necessary is diminished in homo-cellular cultures [95].

The study proved that regulation of the gene expression profile is only a part of the global cancer-induced transcriptomic alteration as illustrated in Figure 5 for 51 RNA polymerase II genes and their binding partners. We believe that changes in the control of genes’ expression fluctuations and their inter-coordination should be also considered and eventually incorporated into the more general measure like the Transcriptomic Distance.

On average, REC increased in T by 2%, in Φ by 28% and in Θ by 20%. An increased control limits the fluctuations of the expression level within a narrow interval, while a decreased control allows the genes to easily adapt their expression levels to various environmental conditions. We believe that the Relative Expression Control is an indirect measure of cell’s priorities to ensure the right expression levels of critical genes for its survival and proliferation. As such, the regulation of REC gives important information of how cancerization has changed the cells’ priorities. We found huge increase of the expression control of POLR2B in both the anaplastic (ΔREC = 485) and the papillary (ΔREC = 206) thyroid cancer cell lines but not in the tumor (ΔREC = 5). The differences might be explained by the potential effects of the encoded protein on the tumor growth as recently reported in glioblastoma [96]. On the other hand, the much stricter control of TAF6 in BPAP cells while looser in the 8505C cells and practically not affected in the surgically removed specimen may indicate different anti-tumor molecular mechanisms [97]. It is interesting to note that POLR2B was up-regulated in both cell lines but not in the tumor (x^(T) = 1.15; x^(Φ) = 74.40; x^(Θ) = 80.06), while TAF6 was down-regulated in both cell lines but not in the tumor (x^(T) = -1.02; x^(Φ) = -1.91; x^(Θ) = -3.24).

For the gene set in Figure 5b, on average, the coordination degree decreased in T (ΔCOORD) = -7.76) but increased in the two cell lines (ΔCOORD) = 8.96 in Φ and ΔCOORD) = 16.80 in Θ). Increased coordination means more synchronization of the genes’ expressions with direct consequences on the pathway efficiency. By contrast, decreased coordination makes the pathway more flexible through letting genes to fluctuate their expression levels more independently. It seems natural to assume that the transcription pathways need more flexibility to adapt to the tissue heterogeneity while in a homo-cellular dish the pathway efficiency comes first.

The high hands out of POLR2L to the transcriptomic alteration in all three types of thyroid cancer samples as measured by the Transcriptomic Distance (Figure 5c) indicate that it might be a good target for anti-thyroid cancer gene therapy as it was suggested for prostate cancer [98]. Interestingly, the main contribution to the Transcriptomic Distance separating the state of a gene in a cancer sample with respect to the normal thyroid may come from any of the three types of individual gene characteristics. Thus, in the case of POLR2L, TD^(T) = 145 comes mainly from the regulation of the expression control (ΔREC^(T) = -142), while in the cell cultures the large TDs (TD^(Φ) = 497; TD^(Θ) = 495) from WIR (WIR^(Φ) = -487; WIR^(Θ) = -452), because of the substantial down-regulations (x^(Φ)= -27.91; x^(Θ)= -25.97). The differences between the tumor and the cell cultures indicate that, although POLR2L is a major player in developing of several forms of thyroid cancers, the involved molecular mechanisms are distinct. Therefore, the therapy targeting this gene should be tailored on the cancer specificity.

Figure 6 and Figure 7 show how cancer remodels the gene network responsible for the DNA Replication and that coupling the genes of the Polymerase II complex to the basal transcription factors in the investigated three types of thyroid malignancy. Expression levels of synergistically expressed genes fluctuate in phase (simultaneously going up and down) across the biological replicas, while those of the antagonistically expressed fluctuate in anti-phase (when one goes up the other goes down). Thus, COR analysis measures indirectly the expression synchronization of the coordinated genes. On the other hand, expression independence tells that most likely there is no direct interaction of the encoded products by the two genes. Missing continuous red/blue or dashed black line between two genes does not mean that they do not interact but that their expression correlation is not statistically significant in that condition. In addition to the substantial yet different remodeling of the transcriptomic network, it is also important to note the very different gene expression profile in the three types of thyroid cancer samples.

There are four main takehome findings in Figure 6:

i) Cancer remodels the DNA replication pathway as can be seen by comparing the gene network in the normal tissue (Figure 6(a)) and the malignant part (Figure 6(b)) of the same thyroid. Changes of the DER genes interaction might provide a complementary explanation of cancer development beyond the driving mutations [99].

ii) Expression correlation is independent of the regulation status of the linked genes. For instance, in “DPD –> Final” there is a significant synergism between the not regulated FEN1 and REC1 in the PTC nodule of the tumor ((COR = 0.978; Figure 6(b) but also between the up-regulated FEN1 and the down-regulated REC1 in the BCPAP cells (COR = 0.951; Figure 6(c)). The FEN1 - REC1 synergism was very near the 0.950 statistical significant cut-off in the normal tissue (COR = 0.944) but less close in the 8505C cells (COR = 0.749). FEN1 was reported to promote cancer cell proliferation, migration and invasion [100] in biliary tumors.

iii) The remodeling depends on the form of the thyroid cancer (see differences between the papillary (Figure 6(c)) and anaplastic (Figure 6(d)) cell lines. For instance, the sygnificant antagonism of the significantly down-regulated genes RPA2 and POLE4 in the BCPAP cells is switched into a significant synergism in 8505C cells (present also in the normal tissue). Although miR-519 was identified as a common upstream regulator of both RPA2 and POLE4 [101], our result might be the first confirmation of their direct interaction.

iv) Each affected person has a distinct DER gene network as illustrated by the differences between the network in Figure 6(b) (PTC collected from the thyroid of a 33y old woman) and that in Figure 6(c) (PTC collected from a 76y old woman [102]). For instance, the independently expressed RPA2 and REC4 in T are syngersically expressed in Φ (and also in Θ). This difference indicates distinct molecular mechanisms of the DNA replication in the papillary cancer nodules of the two PTC donors that requires personalized approach of anti-cancer therapy (e.g. [22,36,48]).

The general findings from Figure 6 can be also retrieved in Figure 7:

i) Cancer remodels the significant expression coordination of polymerases with their targeted partners (compare the correlations in the normal tissue (Figure 7(a)) and the malignant part (Figure 7(b)) of the same thyroid.

ii) Expression correlation is independent of the regulation status of the linked genes. For instance, POLR2H, a prognosis gene for the lung squamus cell carcinoma [101] and ECC3 are synergistically expressed in both T and Φ samples although POLR2H was up-regulated in T but down-regulated in Φ.

iii) The remodeling depends on the form of the thyroid cancer (see differences between the papillary (Figure 7(c)) and anaplastic (Figure 7(d)) cell lines.

iv) Each affected person has distinct gene networks as illustrated by the differences between the correlations in the the two PTC type samples T and Φ. For instance, POLR2H and CDK7 are synergistically expressed in T but antagonistically expressed in Φ.

5. Conclusions

By providing the most theoretically possible comprehensive measure of the transcriptome alteration in cancer, the GFP approach is the alternative of choice to the biomarker paradigm [103]. However, GFP findings should be completed by better equipped laboratories through pathophysiological experiments to identify the molecular mechanisms responsible for controlling the transcripts’ abundances and formation of gene networks.