Two Cohorts, One Network: Consensus Master Regulators Orchestrating Papillary Thyroid Carcinoma

1. Introduction

Thyroid cancer is the most common type of endocrine cancer, making up about 1.0–1.5% of all new cancer diagnoses each year. Its prevalence has steadily increased worldwide, except in Africa, where detection may be limited.[1] At the molecular level, cancer develops due to genetic dysregulation, causing changes in gene expression.[2,3] As high-throughput technologies evolve rapidly, it becomes crucial to analyze and use methods that enable efficient exploration and utilization of genomic data. In this context, systems biology supports integrating diverse information, allowing a comprehensive understanding of the disease on a global scale.

At the molecular level, the progression of thyroid cancer is driven by genetic dysregulation. Some transcription factors (TFs) are described as transcriptional master regulators (TMRs), which control networks of target genes that determine cell differentiation, regulate signaling pathways, modulate hormonal responses, and influence tumor invasiveness, positioning themselves at the top of a regulatory hierarchy.[4]

In this context, Gene Regulatory Network (GRN) inference enables the identification of key TMR-target genes in the development of this disease. GRNs are obtained using computational algorithms that identify the regulatory actions of TFs on all genes in the expression profile, defining the regulon (the set of genes regulated by a TMR) for each TF. This process is known as deconvolution.[5] In this study, we employed network deconvoolution using Mutual Information, which measures dependencies between two variables and was developed based on Shannon’s information theory.[6,7] This type of network infers links between pairs of nodes based on the information they share.[8,9]

In papillary thyroid carcinoma (PTC), there is still a lack of an integrated framework to robustly identify TMRs and regulon networks that control the molecular programs of the disease. Specifically, the hierarchical organization of these regulators and the collective biological functions they influence have not been clearly defined, which limits understanding of tumor progression and reduces the potential to use them as prognostic biomarkers or therapeutic targets[3,4,10] . Our study uses a transcriptomic meta-analysis combining independent cohorts to infer regulatory networks, identify key transcriptional master regulators, and examine their effects on important processes like differentiation, signaling, and tumor plasticity in thyroid cancer.

We primarily focused on identifying transcriptional master regulators (TMRs) and their influence on genes that control the Hallmarks of Cancer. These hallmarks include ten biological traits seen during tumor development.[3,11,12] We identified TMRs that are key in papillary thyroid cancer, including TEAD4. These TMRs not only govern essential transcriptional programs but are also linked to signaling pathways such as NOTCH, MAPK, PI3K/AKT, and estrogen response pathways, as well as processes like epithelial-to-mesenchymal transition and tumor plasticity.

Characterizing these networks and regulatory nodes improves our understanding of thyroid cancer biology, providing potential prognostic biomarkers and therapeutic targets. Zinc-finger enrichment among TMRs. Grouping the 50 TMRs by transcription factor family revealed an overrepresentation of zinc-finger proteins—including ZNF, ZBTB, ZKSCAN, ZFP, and CCCH-type (ZC3H) subfamilies—with eight significant members in our set (ZBTB2, ZBTB25, ZBTB5, ZKSCAN3, ZKSCAN5, ZFP62, ZNF548, ZC3H8). This superfamily is one of the largest in the cell and participates in various cellular processes such as development, differentiation, metabolism, and autophagy; its proteins can interact with DNA, RNA, poly(ADP-ribose), and other proteins.[13] Although zinc-finger factors have been implicated in multiple cancers and several members act as oncogenes[14] , their specific roles in thyroid cancer progression are still not fully understood. Notably, ZBTB2 overexpression can promote proliferation, migration, and invasion.[15] The observed enrichment of zinc-finger TMRs here suggests a potential regulatory hub relevant to papillary thyroid carcinoma progression.

2. Results and Discussion

2.1. Meta-Analysis of Master Regulators

The meta-analysis combining VIPER-based MRA from TCGA-THCA and GSE33630 (papillary thyroid carcinoma, adenocarcinoma, NOS) identified 50 transcriptional master regulators (TMRs) with FDR-adjusted Fisher’s combined p < 0.05, supported across both cohorts. As shown in Figure 1, many transcription factors did not reach strong significance; however, this set of 50 TMRs also showed positive mean NES values (greater than 1), emphasizing their role as key regulators in PTC.

Forkhead (FOX) family TMRs. We identified three TMRs from the Forkhead family. Notably, FOXE1 is a well-characterized regulator of thyroid development and differentiation; its dysregulation has been linked to thyroid tumorigenesis and metastatic progression through modulation of the epithelial–mesenchymal transition via ZEB1.[16,17,18] We also detected FOXP2, which has been implicated in tumor establishment, with overexpression associated with increased proliferation, migration, evasion of apoptosis, metastasis, and immune cell infiltration. In our analysis, FOXQ1 was overexpressed, a pattern reported to correlate with worse prognosis in thyroid cancer patients.

Nuclear receptor TMRs. We identified four transcription factors from the nuclear receptor superfamily—RARA, RARB, RXRG, and ERRγ (also known as ESRRG). This finding is clinically significant because retinoids have been shown to induce redifferentiation of thyroid cells and restore radioiodine uptake, although responses vary; in this context, RARβ/RXRγ expression has been linked to increased treatment sensitivity in thyroid cancer cell lines and tumors.[19,20,21] Additionally, ERRγ (ESRRG) has emerged as a promising drug target, as ERRγ inverse agonists boost sodium–iodide symporter (NIS) expression and improve radioiodine uptake in models of papillary and anaplastic thyroid carcinoma.[22,23,24] This TMR hub and its regulons provide a mechanistic framework for their role in disease, highlighting potential therapeutic opportunities.

ETS family. TMRs include ETV1, ETV4, and ETV5—members of the ETS (E26 transformation-specific) family—who are involved in activating the canonical MAPK pathway in PTC.[25] Their presence among the top 50 TMRs aligns with disease biology. Functionally, ETV4 promotes proliferation and invasion, while ETV5 mediates BRAF^V600E signaling and drives growth and epithelial–mesenchymal transition, highlighting their roles as pro-oncogenic nodes.[26,27,28] Additionally, 15 other transcription factor families are represented in our analysis, each with fewer TMRs; the mapping between the 50 TMRs and their families is summarized in Table 1.

2.2. Putative TMR–TMR Regulatory Interactions

Under the classical definition, a transcriptional master regulator (TMR) is a transcription factor that determines cell fate and typically sits at the top of the regulatory hierarchy.[4] At this level, a TMR can control many downstream genes, while relatively few genes regulate its own activity.[4,29] Additionally, TMRs may regulate the expression of other TMRs, creating transcriptional activation cascades that enhance coordinated responses across gene regulatory networks.[10]

In this context, our goal was to determine whether, within the inferred regulon network, some TMRs not only regulate their downstream targets but also control other TMRs, positioning them as potential key regulators in the transcriptional circuitry of this neoplasm. After identifying the 50 TMRs that were significant and concordant across the two cohorts (TCGA-THCA/GSE33630), we scanned promoter regions of these 50 TMRs on hg38 using the asymmetric TSS window specified in Methods (−2000/+200 on the positive strand; −200/+2000 on the negative strand). Motifs were available for a subset of TMRs based on public databases, and we used those to search for predicted binding sites among the 50 TMRs, thereby assembling a map of putative TMR–TMR connections (Figure 2). We also examined the differential expression patterns for the 50 TMRs and found directional concordance across datasets in all cases (outer ring, blue/red), which reinforces the biological coherence of their behavior. Finally, to place these predictions in the context of thyroid cancer regulation, motif–promoter hits were filtered to retain only those supported by mutual-information edges from our ARACNe-AP deconvolution, providing an integrated view of how these TMRs may coordinate gene expression programs in PTC.

2.3. Functional Relevance of the 50 TMRs

To evaluate the biological importance of the 50 identified TMRs in thyroid cancer, we performed enrichment analysis of the combined regulon network against KEGG pathways and MSigDB Hallmarks (Figure 3). As expected, thyroid hormone synthesis was significantly enriched in the KEGG pathways. We also observed enrichment of valine, leucine, and isoleucine biosynthesis and degradation pathways. Although these amino acids do not directly regulate thyroid function, they are essential for metabolic balance and nutritional health, which can indirectly influence thyroid physiology. Importantly, dysregulation of branched-chain amino acid (BCAA) metabolism has been observed in various cancers.[30] In thyroid cancer, this dysregulation can activate mTORC1 signaling, promoting protein synthesis and cell growth.[31] Metabolomic studies also suggest that BCAAs could be potential biomarkers.[32]

Cytoskeletal remodeling and junctional integrity. Actin cytoskeleton remodeling and tight junctions also appeared as enriched categories among genes regulated by the TMR meta-regulon (Figure 4). This supports our previous finding that several TMR families are involved in epithelial–mesenchymal transition (EMT) and metastasis. Both cytoskeletal reorganization and junctional remodeling become dysregulated in cancer and are essential for executing EMT[33] , as well as for forming invadopodia that promote invasion and metastatic spread.[34]

Estrogen-responsive programs in Hallmarks. Enrichment of MSigDB Hallmarks identified two categories related to estrogen signaling—estrogen response (early) and estrogen response (late) (Figure 4b). It is well established that estrogens act as growth factors in both benign and malignant thyroid cells, activating the MAPK and PI3K pathways, which can coexist with oncogenic mutations such as BRAF mutations or RET/PTC rearrangements.[35] In cell models, estrogen exposure increases VEGF expression, promoting angiogenesis and migration,[36] possibly through an ERα/KRT19 axis.[37] Our findings not only support the importance of these pathways but also expand the network of genes and TMRs that may regulate this program in PTC.

2.4. Pathway Involvement at the Individual TMR Level

After establishing the global functions of the TMRs, we then evaluated their specific roles within signaling pathways through per-regulon ORA for each of the 50 TMRs. Significant enrichment was found for 18 TMRs in KEGG and for 19 TMRs in the MSigDB Hallmark collection. The results of these individual enrichments are summarized in Figure 4.

2.4.1. The Estrogen Response is the Most TMR-Rich Program

Among Hallmark pathways, estrogen response showed the highest number of associated TMRs in both early and late phases (Figure 5). For the early estrogen response, we found significant enrichment of nine TMRs: SMAD9, GRHL3, MAFB, PLAG1, ETV5, TFCP2L1, TEAD4, SETBP1, and EBF4. In contrast, during the late estrogen response, the first five regulators—SMAD9, GRHL3, MAFB, PLAG1, and ETV5—were still present, but TFCP2L1, TEAD4, SETBP1, and EBF4 were not; instead, GATAD2A appeared. This pattern indicates that these five shared TMRs form a core module driving overall estrogen signaling dysregulation, while TFCP2L1, TEAD4, SETBP1, and EBF4 may mainly influence immediate early responses. Consistent with this, early responses in thyroid models involve inducing VEGF, KRT19, and activating MAPK/PI3K signaling, while the late phase includes ongoing angiogenesis, survival pathways, and microenvironmental adaptation—features likely regulated by GATAD2A.

2.4.2. Four TMRs Dominate Pathway Involvement

FOXQ1, RUNX2, SMAD9, and GRHL3, each exhibiting the highest number of enriched categories (four each, Figure 5). In the TMR–TMR promoter analysis, FOXQ1 is regulated by ETV5 and forms a reciprocal loop with GRHL3 (FOXQ1 ↔︎ GRHL3), indicating a feed-forward/feedback module that could position FOXQ1 above or alongside GRHL3 in regulating p53 signaling, apical junctions, and early/late estrogen responses. The RUNX2-centered regulatory node is overexpressed and, in our analysis, controls KRAS signaling, apical junctions, and coagulation—all categories linked to epithelial–mesenchymal transition (EMT), which also appears enriched. Notably, in the ARACNe- and motif-supported network, RUNX2 is regulated by BHLHE40, one of the potential central nodes highlighted by our hub analysis. Our analysis also shows that SMAD9 exhibits significant enrichment for the TGF-β signaling pathway; an intriguing result, given that canonical TGF-β mediators in cancer are SMAD2/3, whereas SMAD1/5/9 usually mediate BMP signaling within the same superfamily.[38,39] This suggests a non-canonical role for SMAD9 in TGF-β signaling in PTC, thereby expanding the repertoire of transcriptional regulators involved in disease progression. Along these lines, the combined enrichment of estrogen response (early/late) and UV response indicates potential crosstalk among hormonal, genotoxic-stress, and TGF-β superfamily pathways.[40,41] Consistent with its central role, SMAD9 is regulated by TEAD4, another potential central node, in the ARACNe/motif-supported network. Overall, these findings highlight the importance of SMAD9 beyond its traditional role in BMP signaling, suggesting it may serve as an alternative regulator in thyroid cancer biology.

2.5. Top-of-Cascade TMRs

We introduce an objective, quantile-based criterion to prioritize “top-of-cascade” candidates on the ARACNe-supported, motif-backed network: transcription factors (TFs) with out-degree in the upper quartile (Q3) and in-degree in the lower quartile (Q1). Under this rule, five TFs emerge as robust driver-node candidates—PBX4 (out/in = 10/3), GATAD2A (7/2), BHLHE40 (7/3), HEY2 (7/3), and TEAD4 (7/3). Their excess of outgoing and paucity of incoming connections suggest a high hierarchical position and a potential capacity to initiate or propagate transcriptional programs in PTC.

This is a conservative assessment because the filter requires both sequence-level motif evidence and ARACNe support: (i) not all panel TMRs have available motifs, and (ii) we used a single representative motif per TF, which does not account for all recognition variants of each regulator. As a result, some plausible interactions may be excluded, and TMRs with high global out-degree could fall below the quantile thresholds under these constraints. Below, we analyze each factor by considering network placement, cross-cohort differential expression, regulon-level enrichments, and prior evidence, showing how these nodes could initiate or propagate transcriptional programs in PTC.

2.5.1. PBX Homeobox 4 (PBX4)

PBX4, a PBX homeobox transcription factor that acts as an HOX cofactor, has been linked to various tumor types. An extensive pan-cancer study reported widespread dysregulation of PBX4, with overexpression in several cancers but notable underexpression in thyroid carcinoma (THCA), along with associations to stage, methylation, immune infiltration, and pathways connected to MAPK signaling and transcriptional misregulation.[42] In our combined THCA analysis (TCGA-THCA and GSE33630), PBX4 was also downregulated in tumors compared to normal tissue, consistent with the pan-cancer findings. Although our master-regulator meta-analysis identified other transcription factors as primary network drivers, the consistent repression of PBX4 in PTC indicates a context-dependent role that should be monitored, including as a potential indicator of pathway rewiring rather than a primary upstream regulator in this disease.

2.5.2. GATA Zinc Finger Domain Containing 2A (GATAD2A)

The network-level signal we observe for GATAD2A (p66α) aligns with functional data in thyroid models. In anaplastic thyroid cancer cell lines, shRNA knockdown of GATAD2A suppressed proliferation and colony formation, induced G2/M arrest, and markedly increased apoptosis as shown by Annexin V/7-AAD and sub-G1 assays, along with consistent cleavage of caspase-3 and PARP.[43] Moreover, in papillary thyroid cancer (PTC), a circRNA-driven pathway (hsa_circ_0058124/NOTCH3/GATAD2A) promotes proliferation, migration/invasion, and in vivo tumor growth; the axis depends on NOTCH3–GATAD2A modulation and is proposed as a disease driver in PTC.[44] Mechanistically, GATAD2A can also function as a p53 co-activator in breast cancer models by physically interacting with p53, enhancing p53 binding at target promoters (e.g., BAX, NOXA, GADD45A, PAI-1), and suppressing growth and migration. Higher GATAD2A levels are associated with a better prognosis.[45] Although these findings are from breast systems, they support a broader role for GATAD2A as a chromatin-linked modulator of tumor suppressive transcriptional programs, consistent with the context-dependent behavior of NuRD subunits.

2.5.3. Basic Helix-Loop-Helix Family Member e40 (BHLHE40)

In our integrated analysis, BHLHE40/DEC1 is overexpressed and sits upstream of RUNX2, whose regulon converges on KRAS signaling, apical junctions, and coagulation—programs tightly linked to EMT—indicating BHLHE40 as a coordinator of pro-oncogenic modules in PTC. This placement aligns with thyroid cancer experiments, which show that BHLHE40 sustains proliferation and invasiveness and functionally cooperates with NOTCH1. Specifically, BHLHE40 regulates NOTCH1; its overexpression accelerates growth and invasion, while γ-secretase inhibition (DAPT) or NOTCH1 silencing reduces the BHLHE40-driven phenotype.[46] Complementing these data, cohort analyses in THCA report BHLHE40 upregulation associated with immune infiltration (B cells, CD4+ T cells, neutrophils, macrophages, and dendritic cells), and functional evidence indicates that BHLHE40 knockdown decreases proliferation, migration, invasion, and metastasis traits.[47] Mechanistically, BHLHE40 can act as either a repressor or an activator, depending on its partners and the chromatin context. In immune and other systems, it modulates cytokine output, cell-cycle control, and metabolism, providing a framework for coupling stress or inflammatory cues to tumor-intrinsic programs in PTC.[48]

In anaplastic thyroid carcinoma (ATC), an H19–miR-454-3p–BHLHE40 axis activates PI3K/AKT signaling: H19 sponges miR-454-3p to de-repress BHLHE40, and BHLHE40 overexpression rescues the anti-proliferative and anti-migratory effects of H19 knockdown.[49] Along with the DEC1↔NOTCH1 positive feedback[46] , these findings outline two converging routes—NOTCH and H19–miR-454-3p–BHLHE40 (PI3K/AKT)—through which BHLHE40 can connect with EMT-related pathways, consistent with our RUNX2-focused pathway observations.

Beyond thyroid cancer, BHLHE40 is widely involved in various malignancies, highlighting its importance in oncology: breast,[50,51,52] colorectal,[53,54,55] choriocarcinoma,[56] endometrial,[57,58] esophageal,[59] gastric,[60,61,62] glioma,[63,64,65] HeLa models,[66,67] hepatocellular carcinoma,[68,69] non–small-cell lung,[70] oral squamous cell carcinoma,[71,72] and pancreatic cancer.[73,74,75,76,77]

2.5.4. Hes Related Family bHLH Transcription Factor with YRPW Motif 2 (HEY2)

In our integrative analysis of papillary thyroid carcinoma (PTC), HEY2 was found to be overexpressed, and its regulon was enriched for estrogen-responsive programs, indicating crosstalk between Notch signaling and hormone-responsive circuitry in PTC. Independent functional data in PTC support a pro-tumorigenic role for HEY2: the miR-599→HEY2 axis directly targets the HEY2 3′UTR (dual-luciferase), and either overexpressing miR-599 or silencing HEY2 reduces proliferation, migration, invasion, and EMT while promoting apoptosis; these effects are associated with increased E-cadherin/Bax and decreased BCL2, N-cadherin, Vimentin, Snail/Slug, and Notch pathway markers (NOTCH1, DLL1, HES1, NICD, JAG1). Overall, these experiments validate HEY2 as a Notch-dependent driver and suggest that restoring miR-599 or inhibiting HEY2 could be a viable therapeutic approach in PTC.[78]

Across various other cancers, HEY2 is often upregulated and linked to aggressive traits, highlighting its pro-oncogenic role. In esophageal squamous cell carcinoma, HEY2 is overexpressed in about 30% of tumors and associates with lymph node metastasis, aligning with Notch pathway activation during disease progression. In non-small-cell lung cancer, a competing endogenous RNA axis involving PRNCR1–miR-448–HEY2 increases HEY2 levels; restoring miR-448 or silencing PRNCR1 reduces proliferation, migration, invasion, and EMT, with partial rescue by HEY2 re-expression. Mechanistically, HEY proteins (including HEY1, HEY2, and HEY-L) are direct Notch targets that act as bHLH transcriptional repressors through E/N-box binding. This framework explains how Notch→HEY2 signaling connects developmental or stress signals to EMT-related outcomes—consistent with the estrogen-response enrichment seen in the HEY2 regulon in PTC.[79,80,81,82]

2.5.5. TEA Domain Transcription Factor 4 (TEAD4)

TEAD4 is a key transcriptional effector of Hippo pathway output, working with YAP/TAZ to regulate programs related to growth, survival, and epithelial plasticity depending on lineage and context. In our integrated PTC analysis, TEAD4 was under-expressed, indicating a relative decrease in canonical YAP–TEAD signaling. Mechanistically, this matches evidence that the tumor suppressor NF2 (Merlin) directly interacts with TEAD4 (FERM and C-terminal regions) to inhibit its palmitoylation, decrease nuclear localization, and promote destabilization—thereby limiting TEAD4-dependent transcription independently of LATS1/2 and YAP.[83] Additionally, an in-vitro thyroid study reports low TEAD4 expression in tumor tissues and cells, demonstrating that forced TEAD4 overexpression reduces viability, migration, and invasion, reverses EMT markers (↓N-cadherin, ↓Vimentin; ↑E-cadherin), and influences Wnt signaling (↑WNT3A). Pharmacologic Wnt blockade (IWR-1-endo) negates these effects, supporting a TEAD4→Wnt axis that limits pro-metastatic programs in thyroid models.[84] A recent review further places TEAD4 at the intersection of epithelial plasticity, stress signaling, and therapy response across cancers, while highlighting lineage-specific behaviors.[85]

Across various cancers, TEAD4 frequently functions as a pro-oncogenic effector of YAP/TAZ; however, its regulatory role and potency are dependent on context, influenced by lineage, cofactors, and Hippo signaling pathways. Notable examples include signals in bladder,[86,87] brain,[88] breast,[89,90,91] cervical,[92] colorectal,[93,94,95,96,97,98] cutaneous melanoma,[99] esophageal,[100,101] gastric,[102,103,104,105,106,107,108] hepatocellular carcinoma,[109,110] laryngeal,[111] lung adenocarcinoma,[112,113] renal cell carcinoma,[114,115,116] and stomach adenocarcinoma.[117] In this context, the PTC-specific downregulation we observe suggests that Hippo pathway regulation—possibly through NF2-mediated TEAD4 inhibition—may counteract other drivers of proliferation in papillary thyroid tumors.[83]

2.6. Therapeutic Implications

Based on our meta-analysis, two key, actionable pathways in PTC emerge: (i) a NOTCH3–GATAD2A module where GATAD2A (p66α) is involved in transcriptional regulation and NuRD-linked chromatin programs, and (ii) a BHLHE40–NOTCH dependency that connects stress and inflammatory signals to tumor-inherent growth and invasion. These findings support biomarker-guided strategies aimed at disrupting oncogenic inputs into NOTCH3–GATAD2A, restoring GATAD2A-driven activation of tumor-suppressing pathways, or targeting vulnerabilities associated with the NuRD complex identified by our TMR subnetwork.[44] Additionally, pharmacologic blockade of NOTCH signaling (such as DAPT) reduces the growth and invasive traits driven by BHLHE40 overexpression, endorsing therapies that target the BHLHE40–NOTCH axis. Since BHLHE40 can act as an activator or repressor by partnering with HDACs, combining modulation of NOTCH signaling with epigenetic regulation presents a promising approach in BHLHE40-high PTC cases.[46,48]

A third druggable node intersects TEAD4 biology. TEAD palmitoylation is crucial for TEAD stability and YAP binding, and small-molecule palmitoylation inhibitors or covalent TEAD ligands can reduce YAP/TEAD transcriptional activity.[87] Additionally, the NF2–TEAD4 interface itself acts as a regulatory chokepoint: NF2 (Merlin) directly binds to TEAD4, inhibiting its palmitoylation, decreasing nuclear residency, and promoting ubiquitin-dependent degradation—thereby limiting TEAD4 activity independently of LATS1/2 and YAP.[83] Given our cohort-specific TEAD4 downregulation, these mechanisms explain why Hippo pathway constraints may be active in PTC and suggest that restoring NF2 function or mimicking its inhibitory contact could suppress TEAD4 where it is oncogenic, while TEAD-axis agents could still target compensatory YAP/TEAD circuitry.

Collectively, these threads support a translational framework that pairs pathway-informed agents with molecular stratification: (a) NOTCH3/GATAD2A-high or BHLHE40-high tumors prioritized for NOTCH inhibition with or without epigenetic co-targeting; (b) YAP/TEAD-signature–high tumors considered for TEAD-directed agents such as verteporfin (YAP–TEAD disruptor), TED-347 (covalent palmitate-pocket binder), or L06 (auto-palmitoylation inhibitor).[85] Prospective testing in PTC models should include functional readouts of EMT, migration/invasion, and Wnt/MAPK outputs, and incorporate response biomarkers for the NOTCH3–GATAD2A, BHLHE40–NOTCH/HDAC, and TEAD/NF2 palmitoylation axes to account for context-specific roles and pathway crosstalk intrinsic to thyroid cancer.

3. Materials and Methods

To ensure robustness and avoid single-cohort bias, we analyzed two independent cohorts of papillary thyroid carcinoma, adenocarcinoma, NOS: the TCGA-THCA project and the GEO study GSE33630.[118] The overall workflow, summarized in Figure 6, combines differential expression analysis, gene set enrichment, regulon inference, and master regulator analysis. The code for this study is available at https://github.com/hachepunto/TiroidesMasterRegulators.

3.1. Data Acquisition and Preprocessing

TCGA-THCA. Transcriptomic data were obtained from the Genomic Data Commons (GDC) using the TCGAbiolinks package. We restricted the cohort to samples with a primary diagnosis of Papillary adenocarcinoma, NOS. From raw count files generated by the STAR–Counts workflow, non-protein-coding features were removed, keeping only genes annotated as protein-coding. Duplicate gene symbols were consolidated by applying the statistical mode across expression values. Genes with at least 10 reads in more than 80% of samples were retained, followed by TMM normalization (trimmed mean of M-values)[119] implemented in edgeR[120], and subsequent transformation to CPM (counts per million).

GSE33630 (GEO). Raw “.CEL” files were downloaded from GEO. Affymetrix Human Genome U133 Plus 2.0 arrays were normalized using RMA (Robust Multi-array Average).[121] Probe annotation was performed with g:Profiler2[122] to map to HGNC symbols, retaining a single identifier per gene by selecting the probe with the highest B-statistic. Genes with low expression (mean log₂ ≤ 4) and low variability (bottom 25% by variance) were subsequently removed.

3.2. Differential Expression Analysis

TCGA-THCA. To identify differentially expressed genes (DEGs) between tumors and normal tissues, we used limma[123] with voom transformation[124] on the CPM-normalized matrix. A linear model was fitted to the tumor versus normal contrast, followed by empirical Bayes variance moderation (eBayes). Genes with adjusted p (FDR) < 0.05 and |logFC| ≥ 1 were considered significant.

GSE33630. Differential expression was analyzed with limma, comparing papillary thyroid carcinoma (PTC) directly to normal samples. After Bayesian moderation, genes were retained based on the same criteria (FDR < 0.05 and |logFC| ≥ 1). Results were consolidated into a single identifier per gene, producing a reliable DEG profile for the GEO cohort.

3.3. Regulon Network Inference

Regulatory networks were inferred using ARACNe-AP (Algorithm for the Reconstruction of Accurate Cellular Networks, Adaptive Partitioning)[125] , based on the curated transcription factor compendium from Lambert et al. (2018).[126]

TCGA-THCA. We used the CPM-normalized expression matrix restricted to protein-coding genes from tumor samples diagnosed as Papillary adenocarcinoma, NOS. The network was generated with a significance threshold of p = 1 × 10⁻⁸, and 100 independent runs with different random seeds were performed, followed by the ARACNe consolidation step to obtain the final network. The resulting network included approximately 234,000 interactions.

GSE33630. Given the smaller sample size (n = 49), we used a more lenient threshold (p = 1 × 10⁻⁴) to account for the reduced statistical power,[6,125] employing the same approach of 100 runs and consolidation. The final network contained about 76,000 interactions, consistent with the cohort size.

3.4. Master Regulator Analysis (MRA)

Master regulator analysis^27,28 was performed using VIPER (Virtual Inference of Protein-Activity by Enriched Regulon Analysis) implementation.²⁹ For each dataset, we used the ARACNe-AP-inferred regulon network along with a tumor-versus-normal differential expression signature. VIPER estimates the regulatory activity of each transcription factor by analyzing the enrichment of its targets within the signature and evaluating significance against permutation-based null models. This approach identifies the transcription factors most likely to drive the transcriptional programs (transcriptional master regulators, TMRs) observed in papillary thyroid carcinoma.

3.5. Integration of MRA Results by Meta-Analysis (Fisher’s Method)

To identify robust and reproducible master regulators across independent cohorts, we integrated the VIPER MRA p-values from TCGA-THCA and GSE33630 using p-value meta-analysis. Each dataset contributed one p-value per transcription factor, derived from the significance of the VIPER enrichment. We combined p-values using Fisher’s method[127] , which transforms them into a log-likelihood statistic:

$${\chi }^2=-2{\sum }_{i=1}^{k}ln\left(p_i\right),$$

(1)

where p_i denotes the p-value for the regulator in dataset i, and y k is the number of studies combined (two in our case). Under the null hypothesis of independence, the resulting statistic follows a χ² distribution with 2k degrees of freedom. We implemented this method using the metap library.[128]

This approach enables the detection of transcription factors that are consistent across both studies, even when the signal within an individual dataset is not sufficiently strong. Additionally, we computed the mean NES and recorded sign concordance across cohorts (without using it as an exclusion criterion). The meta-analysis identified 50 TMRs with FDR < 0.05, which constitute the core of this study.

3.6. Motif Analysis of TMRs

Binding motifs were retrieved from JASPAR2022[129,130,131] as position frequency matrices (PFMs), selecting the latest version for each TF and converting them into position weight matrices (PWMs). We chose JASPAR2022 over newer versions because, at the time of analysis, it was the most recent supported by TFBSTools[132] , ensuring reproducibility and compatibility of the workflow. To expand coverage beyond the initial set of JASPAR-only motifs, we added missing TFs with motifs from CIS-BP³³, standardizing all matrices to PWMs and implementing a source-priority scheme (JASPAR manual > JASPAR automatic > CIS-BP). Metadata for each TF–motif, including source and identifiers, was recorded to ensure reproducibility. For each TMR with an available motif, we scanned the promoters of the 50 TMRs on hg38 using an asymmetric TSS window (−2000/+200 on the positive strand; −200/+2000 on the negative strand) with min.score = 85%. To minimize false positives, we retained only motif–target pairs supported by the ARACNe-AP network in at least one dataset, combining sequence evidence with coexpression/MI support.

3.7. Differential Expression of TMRs Via Meta-Analysis (Fisher + IVW Fixed Effect)

For the 50 validated TMRs, we evaluated whether they show consistent changes in expression between tumors and controls by integrating evidence from both cohorts.

Combining significance levels. For each TMR, differential-expression p-values from each dataset were combined using Fisher’s method.

$${\chi }^2=-2\sum _{i=1}^{k}ln\left(p_i\right),\hspace{1em}\hspace{1em}{\chi }^2\sim {\chi }_{2k}^2,$$

(2)

(k = 2), followed by multiple-testing correction (FDR).

Combining effect sizes (logFC). Complementarily, we combined log fold-changes using a common-effect (fixed-effect) inverse-variance–weighted (IVW) meta-analysis[133,134] :

$${\widehat{\beta }}_{meta}={\displaystyle \frac{{\Sigma }_i{w}_i\widehat{{\beta }_i}}{{\Sigma }_i{w}_i}}, w_i={\displaystyle \frac{1}{{SE}_i^2}}, Z={\displaystyle \frac{{\widehat{\beta }}_{meta}}{\sqrt{\frac{1}{{\Sigma }_i{w}_i}}}},$$

(3)

yielding a two-sided p-value for ${\widehat{\beta }}_{meta}$ that was subsequently FDR-adjusted. Standard errors for each cohort were estimated from limma outputs (t-statistic and logFC).

Finally, we verified directional concordance of the effect (sign of logFC) across datasets. This Fisher + IVW scheme jointly assesses the significance and magnitude of change, improving robustness to cohort-specific variability.

3.8. Functional Enrichment of the Meta-Regulon

To characterize the biological processes collectively controlled by the identified TMRs, we created a meta-regulon defined as the union of targets regulated by the 50 validated TMRs, combining regulons from both datasets. Over-representation analyses (ORA) were performed on this set for MSigDB Hallmarks[135] , KEGG, and GO, for each dataset, using the corresponding universe of expressed genes as background. Term-level p-values were then combined across cohorts using Fisher’s method and adjusted for FDR. This layer summarizes the biological programs shared by the combined network of TMR targets.

4. Conclusions

A cross-cohort meta-analysis of TCGA-THCA and GSE33630 identified 50 transcriptional master regulators (TMRs) in papillary thyroid carcinoma (PTC). These include members of the Zinc Finger, Forkhead, nuclear receptor, and ETS families—groups with established roles in proliferation, migration, metastasis, and thyroid differentiation. Using a quantile-based criterion on the ARACNe+motif network (high out-degree with low in-degree; upper-quartile vs. lower-quartile thresholds), the most consistent top-of-cascade candidates were PBX4, GATAD2A, BHLHE40, HEY2, and TEAD4. Within this framework, we observed the interactions BHLHE40 → RUNX2 and TEAD4 → SMAD9, as well as a reciprocal loop between FOXQ1 and GRHL3, indicating hierarchical modules with potential system-level impact.

Functionally, estrogen response became a key program, with strong Hallmark enrichment observed across cohorts. The early and late components were supported by nine TMRs—SMAD9, GRHL3, MAFB, PLAG1, ETV5, TFCP2L1, TEAD4, SETBP1, and EBF4—highlighting a hormonal axis likely important to PTC initiation and progression.

Regulon-level analyses further implicated SMAD9 in TGF-β signaling, indicating non-canonical activity (usually associated with BMP) and possible crosstalk with estrogen and genotoxic stress pathways. Along with the edges TEAD4 → SMAD9 and BHLHE40 → RUNX2, these findings highlight SMAD9—and the quantile-prioritized TMR set—as key regulators with potential functional and therapeutic significance in PTC.

Limitations. Combining two independent cohorts from different platforms and requiring dual support (motif + ARACNe) enhances robustness but is conservative; it may miss plausible interactions if they fail either filter, and motif coverage (JASPAR/CIS-BP) limits inference. Additional cohorts and functional perturbations will be essential to confirm and extend these findings.This section is not mandatory but can be added to the manuscript if the discussion is unusually long or complex.