This study included eight PTC and four BTN female cases who were treated at The First Affiliated Hospital of Sun Yat-sen University. The diagnosis of PTC was performed according to the WHO criteria. The thyroid specimens and clinicopathologic data were collected after our institutional review board approval.

Primary PTC and BTN were pulverized in liquid nitrogen with mortar. Tissue powder was then cross-linked with 1% formaldehyde and lysed with lysis buffer. Genomic DNA was then sonicated into 200 to 500 bp fragments using ultrasonicator in lysis buffer. The supernatant was diluted, and the lysates were incubated with anti-H3K27ac antibody (ab4729; Abcam, Cambridge, United Kingdom) overnight at 4°C. Protein–DNA complexes were captured with protein A agarose beads. After extensive washing, protein–DNA complexes were eluted and reverse cross-linked. DNA was purified with QIAGEN (Hilden, Germany) PCR purification kit (Cat No. 28106). Ten nanograms of purified DNA was used for downstream library-prep with NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB #E7103; New England Biolabs, Ipswich, MA, United States), following the manufacturer’s instructions. The quality of sequencing libraries was analyzed with bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) and then sequenced by Illumina NovaSeq 6000 platform. Input DNA from each sample were sequenced using different bar codes.

The clean ChIP-seq reads were aligned to the human genome (hg38) using BWA with default parameters. Then, the unmapped reads and non-uniquely mapped reads (mapping quality < 20) were removed by SAMtools, while the PCR duplicate reads were filtered by Picard. The H3K27ac modification regions were defined using MACS. To find out the differential H3K27ac modification region, the modification regions were merged by BEDtools, and the reads of the merged regions were estimated by HTSeq. Finally, DESeq2 was performed to normalize the counts and to detect differential H3K27ac regions (| fold change| ≥ 2 and p value < 0.05).

Primary tissues were pulverized in liquid nitrogen, and RNA was extracted with TRIzol (Invitrogen, Carlsbad, CA, United States) following manufacturer’s protocol. RNA-seq was performed by Berry Genomics Co., China. Briefly, RNA degradation and contamination were monitored on 1% agarose gels, and purity was checked with NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, United States). RNA integrity was assessed with RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). A total amount of 1 μg of RNA per sample was used as input material for the RNA preparations. Sequencing libraries were generated with NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, United States) following manufacturer’s recommendations, and index codes were added to attribute sequences to corresponding sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq 6000 platform.

Differential expression analysis was performed with DESeq; the clean RNA-seq reads were aligned to the human genome (hg38) using HISAT with default parameters. Then, the reads overlapped with genes were estimated by HTseq. The normalization, counts, and detection of differentially expressed genes were performed by DEseq2 (| fold change| ≥ 2 and p value < 0.05).

Genomic Regions Enrichment of Annotations Tool (GREAT)¹ was used to analyze the functional significance of differentiated H3K27ac-modified regions (McLean et al., 2010). Pathway analysis was performed with gene names using DAVID. Selected Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that have p values of less than 0.05 were reported. Pathway enrichment analysis was performed using gene set enrichment analysis (GSEA)² (Liao et al., 2019).

Enhancers and promoters were identified by MACS (Liu, 2014). To identify SEs, all enhancers were ranked according to their total ChIP-seq signal using ROSE (Zhou et al., 2015). Enhancers were sorted and plotted based on H3K27ac signals in ascending order. The cutoff value was set to distinguish SEs from typical enhancers, in which enhancers were assigned as SEs if their H3K27ac signals exceeded this threshold level.

BCPAP and KTC-1 cells were cultured in RPMI-1640 media supplemented with 10% fetal bovine serum (FBS; #10270-106; Gibco, Thermo Fisher Scientific, Waltham, MA, United States). TPC-1 cells were cultured with Dulbecco modified Eagle medium supplemented with 10% FBS (Gibco). All the cells were maintained at 37°C in a 5% CO₂ humidified chamber.

BCPAP, KTC-1, and TPC-1 cells were seeded on day 0 and treated with JQ1 (S7110; Selleck, Houston, TX, United States) or Vorinostat (SAHA) (S1047; Selleck) or vehicle for 48 h with indicated concentrations on day 1. Fresh culture medium with drugs or vehicle was replenished every 24 h.

Total mRNAs were extracted using RNAsimple Total RNA Kit (#DP419; TIANGEN Biotech, Beijing, China). Five hundred nanograms of mRNA was used as template for reverse transcription with PrimeScript RT Master Mix (RR036A; Takara, Dalian, China). cDNAs were then amplified on ABI QuantStudio 5 Real-Time PCR System (Thermo Fisher), and SYBR Green Master Mix (A25742; Thermo Fisher Scientific, Waltham, MA, United States) was used to detect cDNA amplification. GAPDH was used to normalize gene expression. RNA relative expression was calculated using the 2^–ΔΔCT method.

Statistical analyses were performed using R and GraphPad Prism software. Student’s t-test was performed to compare two groups of the continuous variables. All the p values were two-sided, and p < 0.05 was regarded as statistically significant.

To identify histotype-specific landscapes of active chromatin in thyroid tissues, ChIP-seq using H3K27ac (acetylated lysine 27 of histone H3) antibody was performed on eight PTC and four BTN samples from 12 Chinese females (Supplementary Data). Model-based analysis for ChIP-seq (MACS) was used to analyze and identify significant peaks from the samples. Principal component analysis (PCA) indicated that PTC and BTN samples grouped into two clusters according to their genome-wide H3K27ac profiles, suggesting that the differences were etiology specific (Figure 1A). When comparing H3K27ac peaks between PTC and BTN, 395 peaks were unique to PTC samples, while 437 peaks were unique to BTN. With all the differentiated peaks, PTC and BTN samples segregated well on hierarchical clustering (Figure 1B). Differentiated H3K27ac peaks distributed across different regions genome-wide, which showed a greater distribution in intron and intergenic regions with 8.6% PTC-specific and 14.8% BTN-unique peaks located at promoter–transcription start site (TSS) regions (Figure 1C). After the differentiated peaks were annotated with GREAT, we found that gene ontology (GO) terms enriched with PTC-specific peaks included T-cell activation, regulation of lymphocyte activation, leukocyte cell–cell adhesion, and interferon-gamma secretion. Vasculature development, vascular endothelial growth factor signaling pathway, tube morphogenesis, etc., were enriched with BTN-unique peaks (Figure 1D).

Super-enhancers are clusters of enhancers and are critically important in development, differentiation, and oncogenesis (Whyte et al., 2013). To identify SEs, H3K27ac signals from sliding windows containing enhancers were ranked using ROSE (Zhou et al., 2015). We first identified the SEs in the tumor and BTN tissues, and enhancers in each sample were annotated (Supplementary Data). A total of 3,369 SEs were identified in PTC and 3,366 SEs in BTN samples, among which 367 and 364 SEs were specific to PTC and BTN, respectively (Figure 2A). SE genomic sites were categorized into PTC-specific and BTN-specific SE sites. H3K27ac signaling of BTN SEs at PTC-specific genomic sites was significantly lower than that of PTC, contrary to the situation at BTN-specific SE genomic sites (Figure 2B).

Enhancers were ranked based on their H3K27ac signals (Figure 2C). Some SE-adjacent genes in PTC were known to be involved in TC or play a part in oncogenesis, such as EPHB3, ALOX5, and ELF3. In contrast, SEs in BTN included some tumor suppressor genes, which were previously reported in various cancers, such as ARHGAP24 and CA4. SEs were assigned to their putative target genes and were then subjected to pathway enrichment analysis to determine their biological functional classifications (Supplementary Data). As is shown in GO analyses, these SE-adjacent genes of PTC uncovered some immune-related pathways such as immune cell activation and immune response, and other related pathways included positive regulation of ERK1 and ERK2 cascades, and regulation of NF-kappa B transcription factor activity. BTN-specific GO terms revealed pathways of negative regulation of transcription, cell migration, cell growth, etc. (Figure 2D). KEGG pathway analysis indicated some immune-related and cancer-related pathways in SE-adjacent genes of PTC (Figure 2E).

Global transcriptomic features were analyzed in thyroid samples, and histotype-specific patterns of gene expression for the two groups were identified. There were 2,019 differentially expressed genes specific to PTC, while 800 genes showed a higher expression in BTN (Figures 3A,B and Supplementary Data). KEGG enrichment analysis of differentiated genes revealed that genes with higher expressions in PTC were mostly enriched in immune functions such as T-cell receptor signaling pathway, chemokine signaling pathway, antigen processing, and presentation, while BTN-specific genes were over-represented in Rap1 signaling pathway, calcium signaling pathway, mineral absorption, etc. GO annotations also revealed that genes with a higher expression in PTC were mainly involved in pathways such as immune response and inflammatory response. In contrast, BTN samples had a higher expression in metabolism pathways such as oxidation-reduction process, cellular response to hypoxia, and cellular response to zinc ion (Figures 3C,D and Supplementary Data).

The differentiated genes between the two groups were depicted (Figure 3E). Some known oncogenes and tumor suppressor genes were found in PTC-specific or BTN-specific genes, respectively (Figure 3E). GSEA was conducted and showed significant enrichment of PTC-specific genes in cancer-related pathways, such as NF-kappa B, while BTN-specific genes were enriched in metabolism-related pathways (Figure 3F). Besides, GSEA indicated enrichment of immune-related pathways in PTC-specific genes (Figure 3G).

The SEs were assigned to the most adjacent genes, and the expressions of SE-associated genes were higher in the histotype of interest. As exemplified in SE-associated genes in PTC (Figures 4A,C), oncogenes ALOX5 and ELF3 were more highly expressed in PTC compared with BTN tissues. Moreover, The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) data indicated higher levels of ALOX5 and ELF3 in TC compared with adjacent and normal thyroid tissues (Figures 4B,D). The tracks of ChIP-seq and RNA-seq signals of SLCO2A1 and ARHGAP24 were shown to present wide BTN-specific SEs (Figures 4E,G), coinciding with their higher expressions in BTN and in normal thyroid tissues compared with PTC and TC, respectively (Figures 4F,H). Notably, lower transcription level of SLCO2A1 was observed in advanced TC stages and associated with poor disease-free survival (Supplementary Figure 1).

We integrated H3K27ac ChIP-seq and RNA-seq data to map typical enhancers and promoters to their putative targets. As shown in Figure 5A, PTC samples had stronger H3K27ac signals than BTN samples in the genome location of the PLXNC1 gene. After a comparison of active chromatin regions of PLXNC1 with an available database to infer potential enhancers, predicted enhancer–promoter associations were found (Supplementary Figure 2). The RNA-seq revealed that the PLXNC1 was expressed higher in PTC than BTN samples (Figure 5A). Consistently, a higher expression of PLXNC1 in PTC than normal controls was found when analyzing data from TCGA and GTEx samples (Figure 5B). In contrast, CD36 and TRPC6 showed significantly stronger H3K27ac signals and higher expressions in BTN than PTC (Figures 5C–F). Importantly, treatment of thyroid cancer cells with BRD4 inhibitor JQ1 significantly downregulated the expression of PLXNC1 and some PTC-specific genes, while the histone deacetylase inhibitors (HDACi) SAHA treatment induced robust increases of CD36 and TRPC6 expressions in thyroid cancer cell lines (Figure 5G and Supplementary Figure 3). Furthermore, low expression levels of CD36 and TRPC6 were observed in advanced stages of PTC (Figure 5H) and correlated with poor disease-free survival in patients with PTC (Figure 5I), suggesting that these genes are potential predictors for PTC recurrence.