Genetic relationship between Hashimoto`s thyroiditis and papillary thyroid carcinoma with coexisting Hashimoto`s thyroiditis

Hashimoto's thyroiditis (HT) is present in the background of around 30% of papillary thyroid carcinomas (PTCs). The genetic predisposition effect of this autoimmune condition is not thoroughly understood. We analyzed the microarray expression profiles of 13 HT, eight PTCs with (w/) coexisting HT, six PTCs without (w/o) coexisting HT, six micro PTCs (mPTCs), and three normal thyroid (TN) samples. Based on a false discovery rate (FDR)-adjusted p-value ≤ 0.05 and a fold change (FC) > 2, four comparison groups were defined, which were HT vs. TN; PTC w/ HT vs. TN; PTC w/o HT vs. TN; and mPTC vs. TN. A Venn diagram displayed 15 different intersecting and non-intersecting differentially expressed gene (DEG) sets, of which a set of 71 DEGs, shared between the two comparison groups HT vs. TN ∩ PTC w/ HT vs. TN, harbored the relatively largest number of genes related to immune and inflammatory functions; oxidative stress and reactive oxygen species (ROS); DNA damage and DNA repair; cell cycle; and apoptosis. The majority of the 71 DEGs were upregulated and the most upregulated DEGs included a number of immunoglobulin kappa variable genes, and other immune-related genes, e.g., CD86 molecule (CD86), interleukin 2 receptor gamma (IL2RG), and interferon, alpha-inducible protein 6 (IFI6). Upregulated genes preferentially associated with other gene ontologies (GO) were, e.g., STAT1, MMP9, TOP2A, and BRCA2. Biofunctional analysis revealed pathways related to immunogenic functions. Further data analysis focused on the set of non-intersecting 358 DEGs derived from the comparison group of HT vs. TN, and on the set of 950 DEGs from the intersection of all four comparison groups. In conclusion, this study indicates that, besides immune/inflammation-related genes, also genes associated with oxidative stress, ROS, DNA damage, DNA repair, cell cycle, and apoptosis are comparably more deregulated in a data set shared between HT and PTC w/ HT. These findings are compatible with the conception of a genetic sequence where chronic inflammatory response is accompanied by deregulation of genes and biofunctions associated with oncogenic transformation. The generated data set may serve as a source for identifying candidate genes and biomarkers that are practical for clinical application.


Introduction
Hashimoto's thyroiditis (HT), also known as chronic lymphocytic thyroiditis, is induced by an autoimmune response where immune cells produce autoantibodies that target thyroid antigens, in first instance thyroid peroxidase (TPO) and thyroglobulin (TG) leading to atrophy and transformation of thyroid cells [1,2]. This autoimmune reaction is caused by infiltrating autoreactive T lymphocytes and other components of the immune defense system that constitute an inflammatory infiltrate, which under chronic conditions leads to thyroiditis and hypothyroidism [2,3]. Risk factors for HT include environmental factors, such as exposure to ionizing radiation, carcinogenic chemicals, iodine deficiency, and infections, as well as nonenvironmental risk factors, such as genetic susceptibility, that can be linked in certain cases with a history of inherited autoimmune diseases [4,5]. Genome-wide association studies (GWAS) have identified a number of genetic susceptibility loci that are associated with autoimmunity, especially with development of thyroid autoantibodies [6,7]. A systematic review on studies, carried out in their majority on Caucasian populations, reported that incidence of hypothyroidism, which is commonly associated with HT, is 350/100,000/year in women and 80/100,000/year in men [8]. Although controversial discussed, several studies including a meta-analysis on fine-needle aspiration (FNA) and thyroidectomy studies reported HT as a risk factor for papillary thyroid carcinoma (PTC) development, and this observation was mainly based on the proportionally higher incidence of PTC with (w/) HT compared to PTC without (w/o) HT [9, 10].
Thyroid carcinoma (TC) is the most common malignant endocrine cancer with a worldwide increasing incidence and especially PTC as the most common type of thyroid carcinoma contributes to this prevalence [11]. Micro PTCs (mPTCs) are � 10 mm in size and regarded as a low risk and indolent subtype of PTC [12][13][14]. The mortality rates in TC are approximately between 0.20 and 0.60/100,000/year in women and 0.20 and 0.40/100,000/year in men [11]. The female prevalence observed in most of the thyroid lesions is also dominant in PTC w/ HT where cancer incidence is approximately five times higher in females than in males [15]. A meta-analysis on 71 observational studies reported that PTCs w/ HT are more likely to be multifocal but in general have favorable clinicopathological features including less extrathyroid extension, lymph node and distant metastases, and recurrence than PTCs w/o HT [16].
Reactive oxygen species (ROS), which are DNA damaging agents, are a physiological component of thyroid cell function and have elevated expression in HT. For example, a study on public-accessible RNA sequencing (RNA-seq) and microarray data provided genomic evidence that up-regulated genes related to ROS metabolism and apoptosis were significantly more prevalent in PTC w/ HT compared to PTC w/o HT [17]. This led to the notion that higher ROS levels in HT may contribute to PTC development. An immunohistochemistry (IHC) study on PTC-associated tumor markers suggested that in HT, under prolonged exposure to chronic inflammation, a fraction of follicular epithelia show PTC-like nuclear aberrations that might lead to oncogenic transformation [18]. Collaborating with this, a number of studies identified molecular events common to both HT and TC. For instance, loss of heterozygosity (LOH) in the DNA damage repair gene 8-oxoguanine DNA glycosylase (OGG1) has been frequently found in HT and PTCs, but rarely in benign goiter [19]. An IHC study detected, in comparison to normal thyroid tissue, elevated expression of pAKT and the two isoforms AKT1, AKT2 in HT and in well-differentiated TC w/ or w/o HT [20]. In contrast, tumor suppressor PTEN showed positive expression in normal thyroid tissue and HT but only heterogeneous/ absent expression in TC w/ HT. Another IHC study detected positive TP63 expression in HT and PTC, which led to the suggestion that both entities share a pathobiological relationship [21]. Furthermore, ret proto-oncogene (RET)/PTC rearrangements have been detected in a FISH and real-time PCR study at low frequencies in follicular cells of HT, which led to the assumption that these fusions may contribute to early stages of PTC development [22]. In contrast, BRAF mutations that are almost mutually exclusive to RET/PTC rearrangements, are not present in HT w/o coexisting PTC and are comparably more prevalent in PTC w/o HT [23].
The primary aim of the present study was to identify genes and biofunctions, which are shared between HT and PTC w/ HT. Findings should help to better understand the molecular etiopathogenesis between HT and PTC w/ HT and support further assessment of candidate genes for clinical application.

Tumor samples
For microarray expression analysis, 33 thyroid tissue specimens derived from partial or entire thyroidectomies were preserved in RNAlater (Qiagen Inc., Hilden, Germany) and selected for this study according to the histopathology reports. Age of patients ranged between 14 and 67 years with a mean age of 38.7 years, SD = 14.1 (S1 Table). Five of the six micro PTCs (mPTCs) coexisted with HT. The study was approved by the Biomedical Ethics Unit, Faculty of Medicine, King Abdulaziz University, (approval number: 358-10), and the Institutional Review Board of King Faisal Specialist Hospital and Research Center (approval number: IRB2010-07). Written informed consent from the donors or the next of kin was obtained for the use of samples in research. Microarray CEL files from three normal thyroid (TN) samples were downloaded from a website of the microarray vendor (Thermo Fisher Scientific, Waltham, MA) and served as control.

Sample processing
Total RNA was isolated from the samples using the RNeasy Mini-Kit according to the manufacturer's protocol (Qiagen Inc.). RNA concentration was determined using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). RNA integrity was assessed in three cases by gel electrophoretic RNA separation and in 30 cases by using the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). The RNA integrity number (RIN) in these 30 samples ranged between 5.9 and 9.3 (mean, 7.66) and included only RNA in the acceptable quality range between nearly intact and slightly degraded RNA (chosen RIN cutoff value for clinical tumor samples, 5.0). The processed RNA samples were hybridized to HuGene 1.0 ST microarrays according to established protocols [24]. The labeled and washed microarrays were scanned on a GeneChip Scanner 3000 7G and CEL files were generated by the GeneChip Command Console Software AGCC (Thermo Fisher Scientific).

BRAF mutational status
Sequence analysis of BRAF exon 15 on PCR products, by using either a genomic DNA template or a cDNA template, was essentially reported earlier [24,25].

Raw data analysis
The CEL files were imported into Partek Genomic Suite version 6.6 (Partek Inc., Chesterfield, MO) where expression data were normalized using the robust multi-array average (RMA) method in the default settings. A principal component analysis (PCA) was then conducted on the categorized samples and subsequently, lists of significantly differentially expressed genes (DEGs) (false discovery rate (FDR)-adjusted p-value < 0.05 and fold change (FC) > 2) were computed for the four comparison groups HT vs. TN; PTC w/ HT vs. TN; PTC w/o HT vs. TN; and mPTC vs. TN. Statistical assessment of a comparison group HT w/ PTC vs. HT w/o PTC generated no DEGs (FDR-adjusted p-value < 0.05 and FC > 2). From DEGs, represented by multiple probe sets, the set with the highest FC was selected for further processing. To generate intersecting and non-intersecting DEG sets from different comparison groups, a webbased Venn diagram tool was employed (http://bioinformatics.psb.ugent.be/webtools/Venn/) (S2 Table). Where applicable, the Ensemble genome browser and HUGO Gene Nomenclature Committee (HGNC) multi symbol checker were utilized to annotate gene IDs according to the approved symbols [26,27]. An expression dendrogram was generated using the online web tool Heatmapper and employing average linkage clustering and Spearman's rank correlation for distance measurement on the gene set [28]. The generated microarray data set has been deposited at the NCBI's Gene Expression Omnibus (GEO) under accession number GSE138198 and includes renormalized sample files from submission GSE54958.

Pathway and network analysis
Biological significance of DEGs was further interpreted using the Ingenuity Pathway Analysis software and its pathway Knowledge Base (IPA; build version 485516M; Qiagen Inc.). Analysis settings included direct and indirect molecular relationships. Fisher's exact test p-values indicated significant relationships between analyzed data set molecules and frameworks prebuilt or generated de novo by the software application. The Molecule Activity Predictor was employed to predict expression effects of a molecule on other pathway/network molecules as noted in the prediction legends of the respective figures. The overlap percentage in canonical pathways indicates the relative number between uploaded and predefined pathway molecules. Network analysis was performed to explore significance of fit, expressed as a score, between uploaded data set molecules and networks related to specific diseases and functions. The upstream analysis module was employed to assess, by utilizing z-scores, in how far differences in target gene expression are affected by upstream regulators.

Results
To narrow down the genetic relationship between HT and PTC w/ HT, in reference to PTC w/ o HT, mPTC, and TN, we analyzed the molecular expression profiles of a case series consisting of 13 HTs, eight PTCs w/ HT, six PTCs w/o HT, six mPTCs, and three TN samples. Demographic and clinicopathological data of the case series are listed in S1 Table. All tumors, except one, were stage I cases and a V600E BRAF mutation was revealed in eight cases. Similarity of expression profiles of the different groups of cases is presented in a PCA 3D scatter plot demonstrating that expression profiles of the different thyroid disease groups partially overlap whereas TN samples clearly cluster separately (Fig 1). Four comparison groups of DEGs were generated, i.e., HT vs. TN; PTC w/ HT vs. TN; PTC w/o HT vs. TN; and mPTC vs. TN. A Venn diagram displays the number of annotated DEGs, which are either unique for a comparison group or are shared between two or more comparison groups (Fig 2) (S2 Table). In our study, we focused in first instance on the 71 DEGs, which were shared between the two comparison groups HT vs. TN T PTC w/ HT vs. TN. In further instance, we focused on the nonshared 358 DEGs derived from the comparison group HT vs. TN and on the 950 DEGs, which are shared between the four comparison groups, i.e., HT vs.
Top canonical pathways and networks that are most significantly related to the 71, 358, and 950 DEG sets are to varying degrees enriched in immunogenic functions (Table 1).
For all 15 DEG sets, a GO classification in various categories was performed to generate an approximate overview of the functional implications of the sets. The 71 DEGs from the intersection of HT vs. TN \ PTC w/ HT vs. TN harbored the relatively highest number of DEGs related to immune and inflammatory functions; oxidative stress and ROS; DNA damage and DNA repair; cell cycle; and apoptosis ( Fig 3).

Analysis of the 950 DEG set from the intersection of all four comparison groups
In the 950 DEG set, the number up-and downregulated genes was nearly the same. Among genes with highest overexpression were several small nucleolar RNAs, e.g., SNORA71D,

Discussion
This is one of the first microarray expression studies that established gene expression profiles simultaneously in HT, PTCs w/ HT, PTCs w/o HT, mPTCs, and TN samples. By applying a categorization scheme, we were able to generate a comprehensive overview of DEGs shared or non-shared between different types of thyroid lesions with a focus on HT and PTC w/ HT. In the intersection of PTC w/ HT vs. TN \ PTC w/o HT vs. TN we found for instance a strong downregulation of TPO, which is in line with its well-characterized low expression levels in PTC [34]. Although the link between ROS, inflammation and cancer has been extensively discussed, the underlying molecular genetics are not fully explored yet [35]. An expression analysis study on a publicly deposited RNA-seq data set focused on detection of ROS-related genes in PTCs w/ HT vs. PTCs w/o HT and identified 33 genes that were upregulated in PTCs w/ HT compared to nine in PTCs w/o HT [17]. Furthermore, the study found a considerable overlap of the identified upregulated ROS genes with a publicly deposited microarray expression data set. The immune-related expression pattern of genes in PTC was investigated in a number of studies. An protein expression study detected endogenous IgGγ and IgGκ expression in the majority of PTC and colocalization of IgGγ with complement proteins led to the suggestion that immune complexes are formed that may promote PTC cells from escaping the immune response against tumor antigens [36]. The hypothesis that the immunogenic link between HT and PTC may stimulate their synchronous appearance was also discussed in detail [37].  IL2RG is a critical signaling component of a number of interleukin receptors and has been described in a systems biology approach as a factor involved in T cell differentiation signaling during progression of PTC [42]. IFI6 is implicated in apoptosis regulation and exhibits a role in hepatitis C virus infection [43]. The polysialyltransferase ST8SIA4 has been identified in a multiplex meta-analysis of RNA expression datasets as one of the top deregulated genes associated with various immune dysfunctions [44]. STAT1 is implicated in signaling of IFN and a number of interleukins and therefore represents a critical factor of the normal immune response. In addition, STAT1 participates to the morphology and function of the thyroid [45]. An IHC study in HT indicated that tyrosine-phosphorylated STAT1 staining was associated with a thyrocytes to oncocytes transformation as well as with a fraction of infiltrating lymphocytes and germinal macrophages [46]. Of notice, STAT1 constitutes a molecular target for anti-inflammatory treatment [47]. The putative tumor suppressor gene PAG1 encodes a transmembrane adaptor protein that binds to the tyrosine kinase Csk and participates in the negative control of T cell receptor (TCR) signaling [48]. The cytoplasmic enzyme WARS1 catalyzes the ligation of tryptophan to its cognate tRNA during the translation process and beyond this, exerts diverse biological functions in innate immunity, angiogenesis, and IFNG signaling [49]. A microarray mRNA expression profiling study in several major blood cell subtypes derived from normal individuals identified a comparable overexpression of MS4A6A in monocytes/ platelets [50]. The function of MS4A6A in autoimmune-related thyroid lesions remains to be elucidated.
The small GTPase Rac2 is involved in various aspects of host defense including integrin and immunoreceptor signaling, polarization to M2 macrophages and generation of ROS [51][52][53]. Rac2 knockout mice demonstrated a pronounced impairment in tumor growth, angiogenesis, and metastasis [52]. The membrane-bound CYBB is part of the catalytic core of an NADPH oxidase multi-protein complex that produces a phagocytic respiratory burst by generating superoxide in response to pathogens and various other stimuli [53]. The proapoptotic BAK1 protein is a BCL2 family member and an inducer of the intrinsic apoptotic pathway. All three genes, encoding for RAC2, CYBB, and BAK1, were reported to be deregulated among ROS-related genes in an RNA-seq data set analysis between PTC w/ HT and PTC w/o HT [17]. Overexpression of the zinc-dependent metalloproteinase MMP9 has been previously associated with invasion of follicular ML-1 thyroid cancer cells [54]. Furthermore, in vitro assays demonstrated that MMP9 expression is necessary for the serine/threonine kinase ROCK1 to promote PTC cell invasion [55]. TOP2A is a critical factor for sister chromosome segregation and is expressed at highest levels in G2 and M cell cycle phases. The enzyme induces double-stranded DNA breaks, which activate autoimmune regulator (Aire)-associated proteins leading to increased and efficient gene transcription [56]. The serine/threonine protein kinase PLK4 is involved in initiating centriole assembly. Its reversible inhibition may represent a prospective anticancer therapy [57]. Upregulation of BRCA2 is likely a DNA repair response mechanism to elevated DNA damage under chronic ROS exposure. The transcription factor FOXM1 is involved in many crucial biological processes and specifically related to cell cycle progression and cancer evolution [58,59]. FOXM1 overexpression was detected by IHC in 28.4% of PTCs [60]. Further in vitro and in vivo study found that the FOXM1 inhibitor thiostrepton inhibited PTC cell line xenograft tumor growth and this was accompanied by downregulation of FOXM1, MMP9, and MMP2. H2BC14, alias HIST1H2BM, has been identified as one of several molecules that participate in progression through the cell cycle and that are regulated by the p53-p21-DREAM-E2F/CHR pathway [59].

Downregulated genes of the 71 DEG set from the intersection of HT vs. TN \ PTC w/ HT vs. TN
Downregulation of EYA4 expression in cells of epithelial origin is consisting with its tumor suppressor function in these entities [61]. ADGRG7, alias GPR128, is a receptor of the adhesion GPCR family, whose functional relevance in thyroid lesions is virtually unknown; however, promoter studies indicated that ESRα and SP1 play a critical role in positively regulating expression of the receptor [62,63]. ACADL is a peroxisomal enzyme and is part of a family of enzymes known to be regulated by thyroid hormones [64]. Epb41l5 is implicated in the posttranscriptional regulation of cadherin Cdh1 and cell adhesion molecule Itgb1 during epithelial-mesenchymal transition (EMT) and this transition is inhibited by Epb41l5 silencing [65].

Non-intersecting 358 DEG set derived from the comparison group HT vs. TN
The predominant overexpression of IGHV genes, especially IGHV3, that are used for generating autoantibody against TPO is consistent with a review reporting that expression of these genes is associated with HT [66]. Of notice, all IGHV3 genes in our data compilation were included in the 358 DEG set. Similarly, the CD3 components, which were contained in our data compilation, i.e., CD3G, CD3E, and CD3D, were highly overexpressed in the 358 DEG set. In T cells, these molecules generate activation signals when associated in a complex with the TCR. Therapies with CD3 autoantibodies in combination with other therapeutics are envisaged to improve safety and optimize efficacy for treatment of various autoimmune diseases [67]. PDCD1LG2 is expressed from APCs whereas CD28 is expressed from T cells and both immune checkpoint molecules are implicated in inhibitory and stimulatory T cell signaling, respectively [68]. Also of notice, a pronounced number of critical immunomodulators, which were described in an immunogenomic analysis of major TCGA cancer data sets and included PDCD1LG2, BTN3A2, GZMA, BTLA, CD28, ICOS, and IDO1, were contained in the 358 DEG set [69]. The intracellular phosphatase PTPN22 is implicated in immune cell signaling and known as a predisposing gene for autoimmune diseases [70]. A non-synonymous polymorphism in this gene is associated with enhanced susceptibility to a variety of autoimmune diseases [71]. PTPN22 may constitute a target for cancer immunotherapy [70].

The 950 DEG set from the intersection of all four comparison groups
CXCL10 is known as one of the key immunomodulators in cancer and its high expression upon combined INFG and TNFA stimulation has been demonstrated in primary thyrocyte cell cultures [72]. PARP3 is a member of the poly(ADP-ribose) polymerases and is implicated in biological relevant functions including DNA-repair mechanisms and cell cycle regulation [73]. In in vitro experiments, PARP3 expression assisted in TGFB-induced and ROS-supported EMT [74]. PARP3 may gain relevance as a therapeutic target in cancer treatment [75].
An IHC study found a significant association of TP53 positivity in TC with thyroiditis [76]. Increased apoptosis was detected in a case series of HT and thyroid tumors, and it was noted that elevated TP53 and/or CDKN1A protein expression is an indication for possible DNA damage and enhanced apoptosis in HT [77]. Both SOD1 and SOD2, which catalyze the dismutation of superoxide radicals into less reactive hydrogen peroxide and oxygen, were significantly higher expressed in the thyroid lesions compared to TN. A microarray expression analysis of the Oncomine database demonstrated that SOD2 is higher expressed in PTCs and pronounced higher expressed in anaplastic thyroid carcinoma in comparison to normal thyroid tissue [78]. RAD21 is a component of the cohesion core complex that is a critical regulator of sister chromatid cohesion and chromosome segregation. Loss of RAD21 impairs genomic loop domains [79]. Deregulation of RAD21 has been observed in a number of cancers [80]. Different mechanisms of BRCA1 downregulation have been described [81]. Knockdown of BRAC1 in breast cancer cells resulted in increased phospho-Akt and enhanced susceptibility of the cells to PI3K/AKT pathway inhibitors [82]. The secreted prolymphangiogenic factor VEGFD is implicated in tumor growth; however, decreased serum levels of VEGFD were measured in patients with metastatic differentiated thyroid carcinoma [83].

Conclusion
Taken together, this microarray expression study provides a comprehensive overview on DEGs and biofunctions in the thyroid lesions under investigation with a primary focus on those that are shared between HT and PTC w/ HT. Our findings are in line with the conception that in HT a lymphocytic infiltrate is associated with chronic exposure to ROS that promotes DNA damage resulting in oncogenic transformation of thyroid cells. The comprehensive lists of DEGs and biofunctions may serve as a source to assess therapeutic targets and biomarkers, especially in HT and PTC w/ HT.