Different genetic landscapes of papillary thyroid carcinoma and paired benign nodules revealed by integrated multi-omics analysis

doi:10.21203/rs.3.rs-2584167/v1

Download PDF

Research Article

Different genetic landscapes of papillary thyroid carcinoma and paired benign nodules revealed by integrated multi-omics analysis

https://doi.org/10.21203/rs.3.rs-2584167/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Alterations in the genetic landscape of papillary thyroid carcinoma (PTC) compared with coincidental benign thyroid nodules, especially adenomatoid nodules, remain to be demonstrated.

Methods: Multi-omics profiling of whole-exome sequencing, assay for transposase-accessible chromatin using sequencing (ATAC-seq), and transcriptome sequencing were used for analysis.

Results: Chromatin accessibility in the PTC was lower than that in the benign nodules around the transcription start sites (distance <1 kb) with high interpatient heterogeneity of chromatin profiles and distinct open chromatin accessibility. The gene regions around the mutation loci that were only detected in PTC exhibited altered chromatin accessibility between the PTC and benign nodules. Through integrated ATAC-Seq and RNA-Seq analysis, ARHGEF28 and ARHGEF24, genes not previously related to PTC or adenomatoid nodules, were identified. They were overexpressed and hyperaccessible in adenomatoid nodules compared to those in PTC. They were regulated by TEAD4, and hyperaccessible binding sites were enriched in differentially accessible regions in benign nodules. In addition, extrachromosomal circular DNA (eccDNA) analysis derived from ATAC-sequencing showed indolent character, but no PTC-diver genes in the eccDNA was found.

Conclusions: This compendium of multi-omics data provides valuable insights and a resource for understanding the landscape of open chromatin features and regulatory networks in PTC and benign nodule pathogeneses.

Chromatin accessibility

cancer open chromatin landscape

benign thyroid nodule

papillary thyroid carcinoma (PTC)

Thyroid cancer (TC) is the most common endocrine cancer, with an increasing incidence in recent decades [1, 2], and papillary thyroid carcinoma (PTC) accounts for approximately 80–85% of all thyroid cancers [3, 4]. Meanwhile, a dramatic increase in the prevalence of benign thyroid nodules, especially adenomatoid nodules, one of the most common types of hyperplasia, has also been reported [5, 6]. In addition, the coexistence of benign, malignant, and multiple malignant thyroid nodules (TC) has also been observed clinically [7, 8].

TC occurs in approximately 10% of thyroid nodules[7] and can be detected in > 60% of the general population [5]. Only 5% of cases have a clear malignant diagnosis [6]. The key factors involved in the divergence of thyroid follicular epithelial cells into two different lesions from the same original tissue type and the tumor risk assessment (possibility) of benign thyroid nodules for triggering malignancy conversion have been studied [9], even though most thyroid cancers are considered indolent and have a favorable prognosis [10, 11].

Genes carrying driver mutations are critical for tumorigenesis [12]. The tumorigenesis-related genetic profiles of most PTCs have been characterized [12–14], and mutations in genes such as BRAF, NRAS, KRAS, and RET mainly affect disease progression and have clinical significance [15–17]. However, their use for distinguishing genetic alterations between benign and malignant thyroid nodules remains controversial [18]. It has been reported that benign nodules possess a unique molecular signature that differs from PTC and have an independent origin [9]; however, further studies are needed.

Epigenetic profiles, such as DNA methylation [19], histone modifications [20], and noncoding RNA [21, 22] have been characterized in thyroid cancer. Chromosomal accessibility is a powerful method for mapping genome-wide open chromatin patterns in multiple tumors [23, 24] and provides valuable insights into the regulatory mechanisms of TC [25]. However, there was no association between chromosomal accessibility and genomic alterations among patients with paired benign/malignant thyroid nodules. Furthermore, phenotypic changes associated with cancer development likely require widespread changes in gene expression programs that drive invasion, migration, dissemination, and colonization. Therefore, the specific regulatory changes driving changes in gene expression may also play important roles in thyroid follicular epithelial cell differentiation.

In this study, chromosomal accessibility was assessed for paired benign (adenomatoid) and malignant thyroid nodules (from PTC), together with genetic alterations and gene expression, to clarify the association between these lesions. Notably, we can explore the effects of genomic mutations and structural variations on open regulatory elements and understand their association with PTC and the adenomatoid transcriptome. The integration of multi-omics analysis revealed novel insights into the gene regulatory mechanism of thyroid follicular epithelial cell reprogramming.

Patients and clinical information

This study was approved by the Institutional Review Board and Human Ethics Committee of the Air Force Military Medical University (Xi’an, China). Patients with papillary thyroid carcinoma (PTC) and benign thyroid nodules were initially diagnosed at the Tangdu Hospital of Air Force Military Medical University (Xi’an, China) from May to December 2020. Notably, all benign nodules analyzed in the present study were adenomatoid nodules. Information, including patient age, sex, ethnicity, pathology, and tumor stage, was collected for these four patients (Table S1). All patients underwent surgical treatment, and none of them underwent neoadjuvant therapy before surgery. PTC Tumors and matched benign thyroid nodules were obtained during surgery. Two pathologists evaluated all samples to determine pathological diagnosis and tumor cellularity. Written informed consent was obtained from all patients.

Study Design

Integrated multi-omics analysis, including whole-genome sequencing (WGS), assay for transposase-accessible chromatin using sequencing (ATAC-seq) profiles, and transcriptome sequencing, was performed on the PTC and benign nodules. WGS was performed on all four paired samples, and ATAC-seq profiles and transcriptome sequencing were performed on three samples. To understand the landscape of open chromatin features and regulatory networks in PTC and benign nodules, PTC-specific mutations were obtained by filtering the paired benign nodule samples. Differences in chromatin accessibility between PTC and benign thyroid nodules were revealed using the ATAC-seq data. Integrated analysis using RNA-seq data obtained from the same patients identified highly correlated gene expression and ATAC-seq peak intensity. The transcription factors (TFs) identified by ATAC-seq were involved in regulatory roles. The association between gene mutation, chromatin accessibility, and gene expression was also revealed in PTC samples and compared to benign nodule samples.

Atac-seq Sequencing

To profile chromatin accessibility, ATAC-seq was performed using a prokaryotic Tn5 transposase system (Nextera DNA Sample Prep Kit, Illumina). To reduce the contamination of mitochondrial DNA in library preparation, nuclear pellets were first isolated from cells before the Transposition Reaction step. Fresh tissues were homogenized for cell preparation and > 5×10⁴ cell nuclei were obtained for chromatin preparation. A Tn5 incubation time of 30 min at 37°C and tagmented DNA were isolated using the DNA Clean & Concentrator-5 kit. Subsequently, 12–15 PCR amplification cycles were conducted. Quality and quantity were identified using Qubit and Agilent 2100 Bioanalyzers, respectively. Paired-end sequencing (2×150bp) was performed on an Illumina HiSeq 2500 platform to yield an average of 50 M reads/sample. ATAC-seq sequencing was performed by Genergy Biotechnology Co. Ltd. (Shanghai, China).

Atac Data Processing

Data were processed by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). Reads containing adapters, more than 10% unknown nucleotides, and low-quality reads containing more than 50% low-quality (Q-value ≤ 20) bases were removed. Raw sequenced reads were aligned to the GRCh38.p13 reference sequence using Bowtie (version 2). Duplicates were removed using Picard (MarkDuplicates). Peaks were defined using MACS2 with parameters “--shift − 75 --extsize 150 --nomodel --call-summits --nolambda–keep-dup all -p 0.01. The region was defined as a peak when log2 (fold enrichment above background) ≥ 3 and − log10 (p-value) ≥ 3. Using the genomic location information and gene annotation information of the peak, peak-related genes were confirmed using ChIPseeker. In addition, the distribution of peaks in different genomic regions (such as promoter, 5'UTR, 3'UTR, exon, intron, downstream, and intergenic) was determined. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted according to the annotated genes. The MEME suite (http://meme-suite.org/) was used to detect motifs. DiffBind was used to analyze peak differences across groups. Differential peaks were determined to be significant when P༜0.05 and log2 (fold enrichment above background) ≥ 1 in the two comparison groups.

Whole Exon Sequencing (Wes)

WES was performed by Genergy Biotechnology Co. Ltd. (Shanghai, China). Genomic DNA was extracted from formalin-fixed and paraffin-embedded thyroid tissues, including PTC and benign nodules, using a QIAamp DNA Mini Kit (51304, Qiagen). The gDNA met the sequencing requirements of purity (optical density (OD)260/280 > 1.8), and the concentration (50 ng/ml) of each sample was fragmented into pieces with a peak size of ~ 200 bp using a biorupter (Diagenode, Belgium), followed by end-repairing, dA-tailing, and adaptor ligation, as outlined in the KAPA Library Preparation Kit protocol (Kapa, Inc.). The libraries were subjected to minimal PCR cycling and quantified using the Qubit 4.0 Fluorometer (Thermo Fisher, Germany). The PCR products were then subjected to exome sequence capture using Twist 36.6M (TWIST Inc.). The captured libraries were analyzed on an Agilent 2100 Bioanalyzer, and the DNA concentrations were measured using a Qubit 4.0 Fluorometer and then sent for sequencing to generate 2×150 bp paired-end reads using a NovaSeq 6000 platform (Illumina, Inc.).

Wes Data Processing

Sequencing reads were aligned to human genome (UCSC hg 19) using BWA 0.5.9. The aligned BAM file was sorted and merged using Samtools 0.1.19, and the duplications were marked and removed using Picard tools. Then, realignment of all insertions and deletions (INDELs) and recalibration of the base quality were done with Mutect2.The calling of variants, including single nucleotide variants (SNVs) and indels, were accomplished using GATKv3.6. The called variants were filtered by hard filtering and a soft filter (Variant Filtering with VQSR) and further annotated using ANNOVAR. Copy-number variations (CNV) were called using EXCAVATOR2 and annotated using ANNOVAR. SNV, CNV, and indel calling was performed separately for each pair of samples.

Transcriptome Sequencing

Whole transcriptome data were generated using RNA-seq. Total RNA was isolated from tissue samples using the standard TRIzol protocol (Invitrogen, Carlsbad, CA, USA). The RNA quality was examined by gel electrophoresis using a NanoDrop spectrophotometer (Thermo, Waltham, MA, USA). For RNA sequencing, cDNA was generated using Clontech SMART-seq v4 (634888) and libraries were generated using the TruSeq RNA sample preparation kit (Illumina, San Diego, CA, USA). Sequencing was performed using the Illumina HiSeq X Ten instrument by the commercial service of Genergy Biotechnology Co., Ltd. Output sequencer files for each sample were handled by Perl and subjected to a quality control assessment using FastQC v0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The read length was 2150 bp.

Transcriptome Sequencing Data Processing

Clean reads were aligned to the human genome (GRCh38) using STAR with the default settings. The transcript numbers were calculated based on the number of fragments per kilobase of transcript per million fragments mapped (FPKM); the values were calculated using the Cuffnorm program and log2 transformed. Differential expression analysis was performed using the DESeq2 R package, modeling the data with a negative binomial distribution and using Empirical Bayes shrinkage for dispersion and fold change estimation. Functional enrichment analyses were performed using KEGG and GO analyses. P was calculated. All RNA-seq tests were BH-adjusted for multiple testing corrections. Significant genes were categorized based on their functions.

The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) data analysis.

TC RNA-seq normalized count data were analyzed using the Gene Expression Profiling Interactive Analysis (GEPIA) database, including TCGA and Genotype-Tissue Expression (GTEx) data. The GSE64912 and GSE152515 datasets were downloaded from the GEO database. The expression data were log2(x + 1) transformed and quantile-normalized. Differences in TEAD family and target gene expression between TC and normal samples were tested using a two-tailed T-test. Functional enrichment of correlated genes was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) GO Pathway gene sets.

1. The characteristics of gene mutations and chromatin accessibility in co-existing paired benign nodules (PAs) and PTC (CAs)

The CAs were characterized by comparison with PAs, including SNV, indel, and CNV, using WES in four paired benign and malignant thyroid nodules from four patients (Fig. 1A, Table S1). Most SNVs and indels were located in exons and introns, whereas only a small proportion of SNVs were located in intergenic regions (Fig. 1B). However, a total of 718 mutations targeting 657 genes were identified, and the mutation rate (density) differed among the four patients with PTC (Fig. 1C, D; Fig. S1A, B). The transition to transversion ratio of SNV ranged from 0.26–0.73 in PTC (Fig. 1E), while the proportion of variation type T > G and A > C was more correlated with PTC (Fig. S1E). In accordance with a previous report, we identified similar SNVs in the driver genes of PTC, including TTN, BRAF, and MUC16 (Fig. 1F). Mutation signature analysis was performed, and seven enriched mutational signatures (S1, S3, S6, S9, S15, S25, and S28) were identified in PTC samples, in which signatures S1 and S9 were found in all PTC samples (Fig. 1G; Fig. S1F). In terms of CNV, we also observed amplification of known PTC-related genes, including TRDN, COL11A2, and BAGE. Significant deletions of other known cancer-related genes were also identified, such as PRKN. However, no consistent changes in CNV in all of the four samples were found (Table S2, Fig. S1G).

Using ATAC-seq sequencing, we generated chromatin accessibility maps of PTC and coincidental benign nodules (adenomatoid nodules) from three patients with PTC (Fig. 1A). We obtained an average of 97.93% mapability and 61.02 million qualified fragments per sample (Table S3). All samples were enriched for reads at transcription start sites (TSSs) and exhibited the expected periodicity of the insert length (Fig. S2A–C). A global epigenetic landscape of PTC and paired adenomatoid nodules was present (Fig. 2). Chromatin accessibility in CA was lower than that in PA around TSSs (distance < 1 kb), although a broadly similar chromatin landscape was observed (Fig. 2B, Fig. S2D). The number of ATAC-seq peaks (open chromatin regions) in each sample ranged from 36938 to 77749. In total, 45048 high-confidence peaks were identified across all CAs and 44443 peaks were identified across all PAs, representing a total of 18329 distinct peaks (Fig. 2B, Table S3). For instance, a peak associated with CD82 is shared by PA and CA, whereas a peak associated with APOH and FOXP3 is specific to CA, and another peak within PTENP1 is specific to PA (Fig. 2C). We then calculated the overall similarity of the ATAC-seq profiles among all samples using multidimensional scaling, which showed that the samples were clustered into three groups, each one from one patient (Fig. 2D), which was also revealed by PCA analysis (Fig. S2E). Functional enrichment analysis (FEA) by GO showed that peaks were enriched in the proximity of genes related to cell development and morphogenesis specific to PA and CA, whereas the signature common peaks were enriched in the proximity of genes related to cellular metabolism processes (Fig. 2E). The chromatin accessibility of the gene regions around PTC-specific mutation loci is also a concern. Chromatin accessibility around most PTC-specific mutation positions showed obvious differences between PA and CA (Fig. 2F), indicating the potential impact of the mutation on chromatin accessibility.

2. The distinct higher expression of more negative cancer correlation factors in PA associated with hyper-accessibility, including ARHGAP24 and ARHGEF28

Next, we sought to determine the functional consequences of the differentially accessible regions observed in PTC and benign nodule samples. Differences in chromatin accessibility between PTC and benign thyroid nodules were analyzed to explore their impact on PTC and benign nodules. By comparing the signals for each peak in PA and CA, we observed that 72,689 (77.39%) peaks had reduced chromatin accessibility in CA (Fig. 3A). This analysis indicated that a widespread decrease in chromatin accessibility was associated with PTC progression compared with that of benign nodules. We then assessed the differential accessibility between these two groups; a signature set of 144 peaks was identified, including 46 increased and 98 decreased peaks in the PTC group compared to that of paired benign nodules (Fig. 3B). The differentially open chromatin regions were largely TSS distal (> 10⁴ bp), with relatively few proximal promoter regions exhibiting differential accessibility (Fig. 3C and Fig. S3A). Large genomic domains (on the order of 10⁴ to 10⁶ bp) were more likely to have differentially open regions than others regions (< 10³ bp), with different peaks in these domains (Fig. 3C). More domains enriched for increased accessible peaks in CA were located closer to TSS (0-10² bp) than PA (Fig. 3C). A heatmap of the top 50 signature regions highlighted the changes in chromatin accessibility that distinguished the two groups (Fig. 3D). Signature peaks assigned to the CA group were enriched in terms associated with purine and nucleobase-containing compound metabolic processes, whereas those assigned to the PA group were enriched in terms related to the transmembrane receptor protein tyrosine kinase signaling pathway (Fig. 3E).

Using RNA-seq data obtained from the same patients, gene expression in CA samples and PA were examined, and the top 25 genes in each group were clustered in a heatmap (Fig. 4A). Next, we integrated ATAC-seq and RNA-seq data, focusing on chromatin regions and differentially expressed genes. We identified 14 genes that were both differentially expressed and contained significantly differentially accessible peaks (Fig. 4B). Of these peak/gene combinations, three were more accessible and overexpressed in CA samples (right corner), whereas 11 were more accessible with higher RNA expression levels in PA(left corner). Moreover, in the set of signature peaks, we recapitulated differentially accessible peaks surrounding the specific peak annotated genes such as WDR72 (Fig. 4C). However, differential gene expression containing PTC-specific mutations and altered chromatin accessibility between PA and CA was observed only in CCDC40 and SDK1 (Fig. 4D Data RNA-seq).

The more differentially chromatin accessible regions between PA and CA were assigned to the gene regions of WDR72, LRP1B, PCOLCE2, RYR2, ERO1LB, ARHGAP24, CDON, SERTM1, and ARHGEF28. Andtheir expression was significantly higher in benign nodules compared to that in PTC. We further validated their expression using publicly available expression profiling data of thyroid cancers using the GEPIA2 and GEO (Accession No GSE64912, GSE152515) databases, which confirmed their higher expression in normal thyroid tissues compared to that in PTC samples, including ARHGAP24 and ARHGEF28(Fig. S4A B). In addition, CD82 was found to be highly expressed in PTC samples (Fig. S4C). These data indicate that RNA-seq is representative of the expression of these genes, and that the gene expression pattern in benign nodules is similar to that in normal tissues.

3. Differential motif enrichment analysis within the differentially accessible regions revealed TEADs as the key associated transcription factors

By checking whether transcription factor (TF) binding was affected in PTC and/or benign nodules, we determined the enrichment of transcription factor motifs within the differentially accessible regions between PA and CA. The most highly enriched motif was the binding site for the TEAD family of DNA-binding factors in PA (Fig. 5A), while JUN-FOS factors were enriched in CA. Moreover, TEAD4 and TEAD2 showed significantly increased footprints in different accessibility regions when comparing PA to CA (Fig. 5B). TEAD4 expression was further confirmed by RNA-seq (Fig. 5C) and GEPIA2 and GEO (Accession No GSE64912 GSE152515) databases (Fig. S5). This pattern confirmed that chromatin accessibility of TEAD4 target sites was decreased in CA, suggesting that reduced target site binding by benign nodule-specific TFs plays a critical role in the difference between PTC and adenomatoid pathogenesis.

The target genes of TEAD4 and the related pathways were characterized. TEAD4 target genes were obtained from a combined analysis of the hTFtarget database, Gene Transcription Regulation Database (GTRD) database, and RNA-seq data (genes upregulated in benign nodules) (Fig. 5D). A total of 103 genes were identified as potential targets of TEAD4, including ARHGAP24 and ARHGEF28, the expression of which tended to increase in PA compared to CA. In particular, TEAD4 target genes were enriched in cell differentiation and GRPase-mediated signal transduction (Fig. 5E).

4. The additional analysis of extrachromosomal circular DNA (eccDNA) derived from ATAC also showed potential effects on gene expression in CA and coincidental PA

ATAC-seq data can be used to identify eccDNA carrying known driver genes [26]. In this study, eccDNA derived from chromosomes was identified using the new Circle_finder pipeline. A total of 81 eccDNAs were obtained in this study: 48 in PA and 33 in CA (Table S4). The length distribution of eccDNA is shown in Fig. 6A, and 57 eccDNAs were < 1 kb. The eccDNA was derived from almost all chromosomes (Fig. 6B). By pooling all eccDNA identified in this study, the length distribution of these circles revealed characteristic peaks at approximately 400 bases (Fig. 6C). The higher GC content was relative to the genome average (Fig. 6D) and the enrichment of the eccDNA (or microDNA) was mainly located in the gene-up 2000 bp, introns, and intergenic regions (Fig. 6E). Finally, approximately 79% of the small eccDNAs reported here appear to have used flanking sequences of 2- to 15-base microhomology (Fig. 6F) to promote the ligation that gives rise to the circle.

Furthermore, the eccDNA-derived location in the chromatin differed between PA and CA (Fig. 6G). The top five most abundant eccDNAs (tandem gene duplications) identified in each group are listed in Table S4. eccDNA harboring the GPR157 gene (cpgi) was the most abundant eccDNA in the PA samples, while harboring the PROS1 gene (intergenic) in the CA samples. No driver genes for PTC were found in these eccDNAs. Different gene expressions of UPK3B, CCDC110, and SPACA9 were found in the RNA-seq data between CA and PA samples (Data RNA-seq), and no PTC-specific mutations were observed in these genes. Gene ontology analysis of all the genes carried on the eccDNA/duplication loci was conducted (Fig. 6H). These pathways were shown to be related to the cytokine-mediated signaling pathway in PA and DNA replication in CA samples.

Understanding the intrinsic molecular mechanisms that control thyroid follicular epithelial cell differentiation may uncover opportunities for more effective identification and therapy aimed at inhibiting conversion early in the disease course of benign or malignant thyroid nodules to optimize surgical resection and chemoradiation. Here, we analyzed PTC and coincidental benign (adenomatoid) nodules to identify global changes in the genetic landscape between PTC and benign nodules using integrated multi-omics analysis, including WGS, ATAC-seq profiles, and transcriptome sequencing. High heterogeneity of chromatin accessibility among patients was observed, and an unexpectedly dramatic increase in accessibility that occurred in PTC samples was uncovered. Given the heterogeneity of chromatin accessibility, however, ATAC-seq profiles could, to some extent, serve as molecular markers to identify PTC samples from benign nodules.

Integrative analysis combining WES, ATAC-seq, and RNA-seq data from the same patient potentially delineated the effects of genomic alterations on the gene regulatory network between PTC and benign nodules. A previous study showed that somatic point co-mutations in protein-coding genes are closely associated with high-order spatial chromatin folding [27]. Mutant p53 participates in the acquisition of a lethal phenotype by fine-tuning the chromatin landscape [28]. We observed that several genomic regions near the PTC special mutation position showed different chromatin accessibility between PTC and benign nodules. However, the correlation between gene mutations and chromatin accessibility requires further investigation.

It is rational to hypothesize that higher or broader chromatin peaks could be associated with critical genes as they indicate a more active gene regulatory relationship. Several studies have profiled the open chromatin state and its effects on T-lymphoblastic leukemia [29], non-small cell lung carcinoma [30], pancreatic cancer[28], breast cancer [31], and others. Single-cell chromatin accessibility profiling was also performed on glioblastoma and an invasive cancer stem cell population associated with lower survival was identified [32]. Open chromatin character and identified enhancers located within gene bodies could be predictive of correlated RNA and protein expression in TC [25]. The differential expression of genes associated with the differential ATAC-seq peaks was determined. Genes, such as WDR72 [33], LRP1B [34], WWOX [35], RYR2 [36] and CD82 [37] have been reported to be associated with TC. ARHGAP24 and ARHGEF28 have not previously been reported to be associated with PTC or adenomatoid nodules. They are regulated by TEAD4, a single transcription factor that enacts a program of gene expression enriched specifically in differentially accessible regions in benign nodules. TEAD4 target genes, including ARHGAP24 and ARHGEF28, are involved in GTPase-mediated signaling pathways, and accumulating evidence shows that Rho-GTPase-activating proteins (RhoGAPs) play suppressive roles in cancer cell proliferation and metastasis [38].

eccDNAs are somatically mosaic and contribute to intercellular heterogeneity in normal and tumor cells. More recently, it has been shown that circular DNA not only promotes the expression of oncogenes [39], but also amplifies the regulatory regions associated with genes as eccDNA [40]. In this study, ATAC-derived eccDNA-containing genes exhibited no reported PTC-driver genes or PTC-specific gene mutations. This was not consistent with the reported genes derived from ATAC in TC [26], which contained PTC driver genes, such as NARS. However, GO enrichment analysis revealed similar pathways, such as DNA replication in PTC samples, which further proved that chromatin accessibility (ATAC-seq) potentially contributed to the eccDNA. The enriched cytokine-mediated pathways in PA samples may also indicate the critical roles of the inflammatory response in benign nodule processes.

In summary, we generated multi-omics data from PTC and benign nodule samples and demonstrated the potential of this type of data in elucidating different gene regulatory networks. However, this study had some limitations. First, the exact mechanism of genome mutation, chromatin accessibility, and gene expression need to be clarified further, even though alterations in the genome sequence of regulatory elements affected the gene regulatory network of PTC and benign nodules. Second, the sample size was relatively small and was not sufficient to override inter- or intra-patient heterogeneity as much as possible. Further studies with a larger population or incorporating more types of omics data are needed. These results may improve the characterization of the full genetic landscape and regulatory network of PTC and benign thyroid nodules.

ATAC-seq: assay for transposase-accessible chromatin using sequencing

CA: genomic variation of PTC

CNV: copy-number variations

DAVID: Database for Annotation, Visualization, and Integrated Discovery

eccDNA: extrachromosomal circular DNA

GEO: Gene Expression Omnibus

GEPIA: Gene Expression Profiling Interactive Analysis

GO: gene ontology

BP: biological process

GTEx: Genotype-Tissue Expression

GTRD: Gene Transcription Regulation Database

KEGG: Kyoto Encyclopedia of Genes and Genomes

PA: paired benign nodule

PTC: papillary thyroid carcinoma

SNV: single nucleotide variant

TC: thyroid cancer

TCGA: The Cancer Genome Atlas

TF: transcription factor

TSS: transcription start site

WGS: whole-genome sequencing

WES: whole-exon sequencing

Ethics approval and consent to participate

This study was approved by the Institutional Review Board and Human Ethics Committee of the Air Force Military Medical University (Xi’an, China). Written informed consent was obtained from all patients.

Consent for publication

Written informed consent was obtained from all patients.

Availability of data and materials

The data used or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests

Funding

This study was supported by Air Force Military Medical University for funding 2020ZTB11 and Tangdu Hospital for funding 2019JSYJ014.

Authors’ contributions

Guoqiang Bao and Xianli He conceptualized the overall idea, message, and structure. Lijuan Yuan and Ping Yang wrote the bulk of the manuscript and created the figures. All authors contributed ideas, edits, and additional text specific to individual studies described in the manuscript. The authors read and approved the final manuscript.

Acknowledgements

We acknowledge all staff from the department of general surgery in Tangdu Hospital. We would like to thank all the reviewers and the editors for their useful suggestions.

Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin. 2022;72(1):7–33.
Wang J, Yu F, Shang Y, Ping Z, Liu L. Thyroid cancer: incidence and mortality trends in China, 2005–2015. Endocrine. 2020;68(1):163–73.
Lim H, Devesa SS, Sosa JA, Check D, Kitahara CM. Trends in Thyroid Cancer Incidence and Mortality in the United States, 1974–2013. JAMA. 2017;317(13):1338–48.
Fagin JA, Wells SA Jr. Biologic and Clinical Perspectives on Thyroid Cancer. N Engl J Med. 2016;375(23):2307.
Singh Ospina N, Iniguez-Ariza NM, Castro MR. Thyroid nodules: diagnostic evaluation based on thyroid cancer risk assessment. BMJ. 2020;368:l6670.
Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, Pacini F, Randolph GW, Sawka AM, Schlumberger M, et al. 2015 American Thyroid Association Management Guidelines for Adult Patients with Thyroid Nodules and Differentiated Thyroid Cancer: The American Thyroid Association Guidelines Task Force on Thyroid Nodules and Differentiated Thyroid Cancer. Thyroid. 2016;26(1):1–133.
Wong R, Farrell SG, Grossmann M. Thyroid nodules: diagnosis and management. Med J Aust. 2018;209(2):92–8.
Durante C, Grani G, Lamartina L, Filetti S, Mandel SJ, Cooper DS. The Diagnosis and Management of Thyroid Nodules: A Review. JAMA. 2018;319(9):914–24.
Ye L, Zhou X, Huang F, Wang W, Qi Y, Xu H, Yang S, Shen L, Fei X, Xie J, et al. The genetic landscape of benign thyroid nodules revealed by whole exome and transcriptome sequencing. Nat Commun. 2017;8:15533.
Grogan RH, Kaplan SP, Cao H, Weiss RE, Degroot LJ, Simon CA, Embia OM, Angelos P, Kaplan EL, Schechter RB. A study of recurrence and death from papillary thyroid cancer with 27 years of median follow-up. Surgery. 2013;154(6):1436–46. discussion 1446 – 1437.
Albero A, Lopez JE, Torres A, de la Cruz L, Martin T. Effectiveness of chemotherapy in advanced differentiated thyroid cancer: a systematic review. Endocr Relat Cancer. 2016;23(2):R71–84.
Cancer Genome Atlas Research N. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676–90.
Macerola E, Proietti A, Poma AM, Ugolini C, Torregrossa L, Vignali P, Basolo A, Materazzi G, Elisei R, Santini F et al. Molecular Alterations in Relation to Histopathological Characteristics in a Large Series of Pediatric Papillary Thyroid Carcinoma from a Single Institution. Cancers (Basel) 2021, 13(13).
Haroon Al Rasheed MR, Xu B. Molecular Alterations in Thyroid Carcinoma. Surg Pathol Clin. 2019;12(4):921–30.
Zou M, Baitei EY, Alzahrani AS, BinHumaid FS, Alkhafaji D, Al-Rijjal RA, Meyer BF, Shi Y, Concomitant RAS. RET/PTC, or BRAF mutations in advanced stage of papillary thyroid carcinoma. Thyroid. 2014;24(8):1256–66.
Vuong HG, Altibi AMA, Duong UNP, Hassell L. Prognostic implication of BRAF and TERT promoter mutation combination in papillary thyroid carcinoma-A meta-analysis. Clin Endocrinol (Oxf). 2017;87(5):411–7.
Prete A, Borges de Souza P, Censi S, Muzza M, Nucci N, Sponziello M. Update on Fundamental Mechanisms of Thyroid Cancer. Front Endocrinol (Lausanne). 2020;11:102.
Zhang M, Lin O. Molecular Testing of Thyroid Nodules A Review of Current Available Tests for Fine-Needle Aspiration Specimens. Arch Pathol Lab Med. 2016;140(12):1338–44.
Park JL, Jeon S, Seo EH, Bae DH, Jeong YM, Kim Y, Bae JS, Kim SK, Jung CK, Kim YS. Comprehensive DNA Methylation Profiling Identifies Novel Diagnostic Biomarkers for Thyroid Cancer. Thyroid. 2020;30(2):192–203.
Celano M, Mio C, Sponziello M, Verrienti A, Bulotta S, Durante C, Damante G, Russo D. Targeting post-translational histone modifications for the treatment of non-medullary thyroid cancer. Mol Cell Endocrinol. 2018;469:38–47.
Le F, Li HM, Lv QL, Chen JJ, Lin QX, Ji YL, Yi B. lncRNA ZNF674-AS1 inhibits the migration, invasion and epithelial-mesenchymal transition of thyroid cancer cells by modulating the miR-181a/SOCS4 axis. Mol Cell Endocrinol. 2022;544:111551.
Li Z, Wang H, Deng X, Zhang J, Wang L, Tang W, Youa W, Nian W. Construction and Validation of an Immune-Related lncRNA Prognosis Model for Thyroid Cancer.Comb Chem High Throughput Screen2022.
Corces MR, Granja JM, Shams S, Louie BH, Seoane JA, Zhou W, Silva TC, Groeneveld C, Wong CK, Cho SW et al. The chromatin accessibility landscape of primary human cancers. Science 2018, 362(6413).
Polak P, Karlic R, Koren A, Thurman R, Sandstrom R, Lawrence M, Reynolds A, Rynes E, Vlahovicek K, Stamatoyannopoulos JA, et al. Cell-of-origin chromatin organization shapes the mutational landscape of cancer. Nature. 2015;518(7539):360–4.
Sanghi A, Gruber JJ, Metwally A, Jiang L, Reynolds W, Sunwoo J, Orloff L, Chang HY, Kasowski M, Snyder MP. Chromatin accessibility associates with protein-RNA correlation in human cancer. Nat Commun. 2021;12(1):5732.
Kumar P, Kiran S, Saha S, Su Z, Paulsen T, Chatrath A, Shibata Y, Shibata E, Dutta A. ATAC-seq identifies thousands of extrachromosomal circular DNA in cancer and cell lines. Sci Adv. 2020;6(20):eaba2489.
Shi Y, Su XB, He KY, Wu BH, Zhang BY, Han ZG. Chromatin accessibility contributes to simultaneous mutations of cancer genes. Sci Rep. 2016;6:35270.
Zampieri C, Panatta E, Corbo V, Mauriello A, Melino G, Amelio I. p53 mutations define the chromatin landscape to confer drug tolerance in pancreatic cancer.Mol Oncol2021.
Erarslan-Uysal B, Kunz JB, Rausch T, Richter-Pechanska P, van Belzen IA, Frismantas V, Bornhauser B, Ordonez-Rueada D, Paulsen M, Benes V, et al. Chromatin accessibility landscape of pediatric T-lymphoblastic leukemia and human T-cell precursors. EMBO Mol Med. 2020;12(9):e12104.
Wang Z, Tu K, Xia L, Luo K, Luo W, Tang J, Lu K, Hu X, He Y, Qiao W, et al. The Open Chromatin Landscape of Non-Small Cell Lung Carcinoma. Cancer Res. 2019;79(19):4840–54.
Achinger-Kawecka J, Valdes-Mora F, Luu PL, Giles KA, Caldon CE, Qu W, Nair S, Soto S, Locke WJ, Yeo-Teh NS, et al. Epigenetic reprogramming at estrogen-receptor binding sites alters 3D chromatin landscape in endocrine-resistant breast cancer. Nat Commun. 2020;11(1):320.
Guilhamon P, Chesnelong C, Kushida MM, Nikolic A, Singhal D, MacLeod G, Madani Tonekaboni SA, Cavalli FM, Arlidge C, Rajakulendran N et al. Single-cell chromatin accessibility profiling of glioblastoma identifies an invasive cancer stem cell population associated with lower survival.Elife2021,10.
Schulten HJ, Hussein D, Al-Adwani F, Karim S, Al-Maghrabi J, Al-Sharif M, Jamal A, Bakhashab S, Weaver J, Al-Ghamdi F, et al. Microarray expression profiling identifies genes, including cytokines, and biofunctions, as diapedesis, associated with a brain metastasis from a papillary thyroid carcinoma. Am J Cancer Res. 2016;6(10):2140–61.
Prazeres H, Torres J, Rodrigues F, Pinto M, Pastoriza MC, Gomes D, Cameselle-Teijeiro J, Vidal A, Martins TC, Sobrinho-Simoes M, et al. Chromosomal, epigenetic and microRNA-mediated inactivation of LRP1B, a modulator of the extracellular environment of thyroid cancer cells. Oncogene. 2011;30(11):1302–17.
Dias EP, Pimenta FJ, Sarquis MS, Dias Filho MA, Aldaz CM, Fujii JB, Gomez RS, De Marco L. Association between decreased WWOX protein expression and thyroid cancer development. Thyroid. 2007;17(11):1055–9.
Xu N, Zhang D, Chen J, He G, Gao L. Low expression of ryanodine receptor 2 is associated with poor prognosis in thyroid carcinoma. Oncol Lett. 2019;18(4):3605–12.
Chen Z, Mustafa T, Trojanowicz B, Brauckhoff M, Gimm O, Schmutzler C, Kohrle J, Holzhausen HJ, Kehlen A, Klonisch T, et al. CD82, and CD63 in thyroid cancer. Int J Mol Med. 2004;14(4):517–27.
Yang W, Wang B, Yu Q, Liu T, Li T, Tian T, Jin A, Ding L, Chen W, Wang H, et al. ARHGAP24 represses beta-catenin transactivation-induced invasiveness in hepatocellular carcinoma mainly by acting as a GTPase-independent scaffold. Theranostics. 2022;12(14):6189–206.
Wu S, Turner KM, Nguyen N, Raviram R, Erb M, Santini J, Luebeck J, Rajkumar U, Diao Y, Li B, et al. Circular ecDNA promotes accessible chromatin and high oncogene expression. Nature. 2019;575(7784):699–703.
Morton AR, Dogan-Artun N, Faber ZJ, MacLeod G, Bartels CF, Piazza MS, Allan KC, Mack SC, Wang X, Gimple RC, et al. Functional Enhancers Shape Extrachromosomal Oncogene Amplifications. Cell. 2019;179(6):1330–1341e1313.

Download PDF

Version 1

posted

You are reading this latest preprint version

Different genetic landscapes of papillary thyroid carcinoma and paired benign nodules revealed by integrated multi-omics analysis

Status:

Version 1

Abstract

Figures

Background

Materials And Methods

Patients and clinical information

Study Design

Atac-seq Sequencing

Atac Data Processing

Whole Exon Sequencing (Wes)

Wes Data Processing

Transcriptome Sequencing

Transcriptome Sequencing Data Processing

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1