Transcriptome, Methylome and Genomic Variations Analysis of Ectopic Thyroid Glands

Background Congenital hypothyroidism from thyroid dysgenesis (CHTD) is predominantly a sporadic disease characterized by defects in the differentiation, migration or growth of thyroid tissue. Of these defects, incomplete migration resulting in ectopic thyroid tissue is the most common (up to 80%). Germinal mutations in the thyroid-related transcription factors NKX2.1, FOXE1, PAX-8, and NKX2.5 have been identified in only 3% of patients with sporadic CHTD. Moreover, a survey of monozygotic twins yielded a discordance rate of 92%, suggesting that somatic events, genetic or epigenetic, probably play an important role in the etiology of CHTD. Methodology/Principal Findings To assess the role of somatic genetic or epigenetic processes in CHTD, we analyzed gene expression, genome-wide methylation, and structural genome variations in normal versus ectopic thyroid tissue. In total, 1011 genes were more than two-fold induced or repressed. Expression array was validated by quantitative real-time RT-PCR for 100 genes. After correction for differences in thyroid activation state, 19 genes were exclusively associated with thyroid ectopy, among which genes involved in embryonic development (e.g. TXNIP) and in the Wnt pathway (e.g. SFRP2 and FRZB) were observed. None of the thyroid related transcription factors (FOXE1, HHEX, NKX2.1, NKX2.5) showed decreased expression, whereas PAX8 expression was associated with thyroid activation state. Finally, the expression profile was independent of promoter and CpG island methylation and of structural genome variations. Conclusions/Significance This is the first integrative molecular analysis of ectopic thyroid tissue. Ectopic thyroids show a differential gene expression compared to that of normal thyroids, although molecular basis could not be defined. Replication of this pilot study on a larger cohort could lead to unraveling the elusive cause of defective thyroid migration during embryogenesis.


Introduction
Permanent primary hypothyroidism is the most common congenital endocrine disorder. In up to 85% of cases, it results from thyroid dysgenesis, a condition comprised of defects in the differentiation, migration or growth of thyroid tissue. Of these defects, incomplete migration resulting in ectopic thyroid tissue (sub-lingual thyroid) is the most common (up to 80%). The etiological diagnosis is established through thyroid scintigraphy [1]. Ectopic thyroids are smaller (i.e. they lack the lateral lobes that are characteristic of orthotopic thryroids) but are otherwise normal (i.e. they have a normal follicular architecture and their capacity to trap and organify iodine and to produce thyroid hormones and thyroglobulin is intrinsically normal [2,3,4,5]). This suggests that the hypothyroidism of subjects with thyroid ectopy is due to a smaller amount of tissue (hypoplasia), which is a consequence of the migration defect, and not to defects in differentiation or in histological organization of the thyroid follicular cells.
Congenital hypothyroidism from thyroid dysgenesis (CHTD) is a heterogeneous disease, which exists in familial (2%) and nonfamilial (sporadic, 98%) forms [6]. Moreover, the results of a survey of monozygotic twins yielded a discordance rate of 92% [7], which together with the female predominance in CHTD [8] suggest that complex non-Mendelian mechanisms underlie this condition. On the other hand, environmental causes operating in utero are unlikely because: (i) no temporal or seasonal trends for CHTD have been observed [9] and (ii) MZ twins who are discordant for CHTD have similar birth weight (G. Van Vliet, personal communication). Germinal mutations in thyroid related transcription factors NKX2.1, FOXE1, PAX-8, and NKX2.5 have been identified in only 3% of patients with sporadic CHTD and linkage analysis excluded these genes in some multiplex families with CHTD, which is consistent with a complex genetic contribution [10]. Together these findings indicate the involvement of novel genes and pathways and underlines the importance of somatic epigenetic or genetic events [11].
Combining data of gene expression, DNA methylation and DNA copy number has led to the identification of novel genetic regulators of cancer [12,13]. Consistent with this approach, we aim to assess whether the transcriptome of ectopic thyroids is shaped by somatic genomic or epigenomic variations (Figure 1).

Results
Expression array identified 1011 genes that are more than two-fold induced or repressed We used microarray analysis to compare the genome-wide RNA expression profile of normal (orthotopic; n = 1) versus abnormal (ectopic; n = 3) thyroid tissue. We identified 1833 differentially expressed genes, and a total of 1011 genes were induced (n = 522) or repressed (n = 489) more than two-fold. The forty genes with the highest differential expression are listed in Table 1. To validate the differential expression identified by microarray analysis, we performed quantitative real-time PCR (qRTPCR) of 100 genes in ectopic thyroids compared with the same commercial control (Ambion) used for the arrays; these 100 genes included highly differentially expressed genes and genes known to play a role in the thyroid function. Overall, there was a highly significant correlation between microarray and qRTPCR (Pearson correlation coefficient of 0,86, p,2.2 e-16)( Figure S1).
Functional annotation clustering of the 1011 differentially expressed genes showed enrichment for developmental processes To assess whether the differentially expressed genes are related to development and organogenesis, we classified the 1011 differentially expressed genes into gene ontology (GO) groups using DAVID (Database for Annotation, Visualization and Integrated Discovery) with medium classification stringency. Table  S1 shows the top three clusters for the 1011 differentially expressed genes (more than two-fold induced or repressed) with enrichment scores greater than 6 (i.e. p,1 E -06). Two of them are clusters of genes enriched for developmental processes. We next clusterized separately the induced (n = 522) and the repressed (n = 489) genes using DAVID according to GO terms with high classification stringency. The top five clusters of induced genes with enrichment   (Table  S2). Finally, the analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with DAVID shows enrichment (p,0.05) for eight pathways, among which three are associated with cell-to-cell interaction (i.e. MAPK signaling pathway, focal adhesion and cell communication) (Figure 2).
Additional validation against age-matched hyperfunctioning thyroid tissue selects 19 genes whose expression is related to thyroid ectopy and independent of thyroid activation state Level of activity (i.e. thyroid hormone production) and consequently expression of genes involved in thyroid hormone production is sensitive to thyroid stimulating hormone (TSH) [14,15]. CHTD patients have high TSH at diagnosis, after which the level of their TSH may vary depending on the compliance to the treatment; or the TSH is high for years in cases of delayed treatment, as in our case 1 (Table S3). To correct for differential TSH-dependent activation of thyroid tissue, we extended the qRTPCR analysis for the 100 validated genes to three hyperfunctioning thyroid nodules (i.e. somatic activation of the TSH receptor) that had arisen in otherwise normal (orthotopic) thyroids. Of the 100 genes, 19 showed consistently induction or repression when compared to all types of controls (i.e. hyperfunctioning thyroid nodules and the commercial control) ( Table 1). The 81 remaining genes showed divergent expression: i.e. they were either induced in ectopic thyroid when compared to the commercial control but repressed when compared to the three hyperfunctioning thyroid nodules or vice versa (Table S4), suggesting that the expression of those 81 genes was dependent of thyroid activation state.
Pathway analysis of the 19 selected genes showed association with Wnt signalling pathway whereas the 81 excluded genes were enriched for thyroid hormonogenesis and function We then asked whether the use of two different types of control was able to exclude genes that are dependent on thyroid activation state and to identify genes associated exclusively with defective migration of the thyroid. To consider the level of expression (which is not possible with DAVID), we used the Ingenuity Pathways Analysis software for GO annotation and pathway analysis (IngenuityH Systems, www.ingenuity.com). The 19 selected genes (Table 2) were enriched for the Wnt pathway, dendritic cell maturation pathway, and embryonic developmental functions (for each enrichment scores greater than 1.5; p,0.02 and Benjamini-Hochberg multiple correction p,0.05). As expected, the 81 genes (Table S4) excluded because of divergent expression are enriched for thyroid hormonogenesis and function: thyroglobulin (TG), thyroid peroxidase (TPO), deiodinase type II (DIO2), deiodinase type I (DIO1), dual oxidase 2 (DUOX2), paired box gene 8 (PAX8), thyroid stimulating hormone receptor (TSHR) and thyroid hormone responsive SPOT14 homolog (THRSP) (i.e for GO term endocrine system disorders, the enrichment score is 2; p,0.01 and Benjamini-Hochberg multiple correction p,0.01). To control for tissue quality, we also analyzed 10 unregulated genes which are either well described thyroid-related transcription factors (FOXE1, HHEX, NKX2.1, NKX2.5) [16] or genes involved in the Wnt pathway (CTNNB1, GSK3B, CDH1, APC, AXIN1, AXIN2) [17]. CTNNB1 and CDH1 showed divergent expression; APC, AXIN1, AXIN2, FOXE1 and NKX2.1 showed increases of convergent expression; and NKX2.5 expression was massively increased when compared to orthotopic thyroids but only slightly increased when compared to the commercial control, which therefore might be considered as divergent. None of the thyroid related transcription factors (FOXE1, HHEX, NKX2.1, NKX2.5) showed decreased expression regardless of the control used (Table 3).

Differential gene expression in ectopic thyroid is independent of methylation
The next step was to assess whether somatic changes in DNA methylation play a role in dysregulation of gene expression in ectopic thyroids. Methylation profiling by methylated DNA immunoprecipitation (MeDIP) and MeDIP-chip was performed by hybridizing pairs of enriched methylated fraction (IP) and normal fraction (IN) of genomic DNA from our three ectopic thyroids and three controls orthotopic hyperfunctioning thyroid nodules. The methylation profile was similar between the ectopic and orthotopic thyroids: after multiple test correction, there was no statistically significant difference (i.e. no region with a less stringent False Discovery Rate threshold of 0.1) (data not shown). Consequently, no correlation was found between the differential expression in ectopic thyroids and the global methylation profile.

Differential gene expression in ectopic thyroid is independent of thyroid-specific CNVs
To assess whether thyroid-specific (i.e. absent in matched leucocytes) CNVs shape gene expression in ectopic thyroids, we used array comparative genomic hybridization (aCGH) of the ectopic thyroid DNA with matched leucocytes. By analyzing data as described in Materials and Methods, we found four thyroid-specific CNVs (three validated by qPCR), which are reported variants as reported in the Database for Genomic Variants (http://projects. tcga.ca/variation) (Table S5). No correlation was found between thyroid-specific CNVs and differentially expressed genes in ectopic thyroids.

Discussion
Generally, CHTD is sporadic and shows discordance between MZ twins [7]. Somatic genetic or epigenetic events might therefore have a role in the etiology of this condition. The objective of this study was to assess whether somatic molecular events account for the failure of migration of ectopic thyroids. Therefore, we conducted the first integrative analysis of transcriptome, DNA methylation and structural variants (CNV) in ectopic thyroids.
We found altered expression in genes and pathways that might play a significant role in abnormal thyroid development (e.g. Wnt signaling pathway). Interestingly, none of the thyroid related transcription factors (FOXE1, HHEX, NKX2.1, NKX2.5) showed decreased expression, whereas PAX8 expression was associated with thyroid activation state. This is a direct indication that the expression of these known candidate genes is at least neutral in ectopic thyroid and is consistent with the observation that the coding sequences for FOXE1, NKX2.1 and PAX8 were normal in case #1 [5].
Four pathways identified by analysing the results of expression arrays (i.e. focal adhesion, antigen processing and presentation, cell communication, cell adhesion molecules and Type I diabetes) have been identified independently in hyperfunctioning thyroid nodules [18]. However, our results identify mostly induced genes in these pathways ( Figure 1) whereas repression of these genes was observed in the aforementioned study [18]. To obtain, in ectopic thyroid, an opposite expression pattern when compared to that of orthotopic hyperfunctional thyroid (i.e. thyroid with somatic TSH receptor activating mutation) is plausible, but it underlines also the need to consider the differential activation of the TSH-receptor signaling pathways in our samples. Consequently, we have excluded 81 validated genes for which expression was associated with TSH-driven thyroidal activity.
The 19 selected genes whose expression was dependent on thyroid location (i.e. ectopy) were enriched for pathways involved in cellular movement (i.e. Wnt pathway and dendritic cell maturation pathway). This association has biological plausibility especially for the Wnt pathway. First, non-canonical Wnt pathway is crucial for cell migration [19] and development of organs of endodermal origin (e.g. intestine, lung, pancreas) [20]. There is indirect evidence for the involvement of the non-canonical Wnt pathway in the developing thyroid in mice [21], even though the canonical Wnt/beta-catenin pathway seems to be inactive during thyroid development in mice and humans [21,22]. Second, as Wnt signaling is implicated in development and cancer [23], to find an association between Wnt pathway and thyroid ectopy (i.e. failure of proper thyroid migration during development) makes biologically sense. Indeed, SFRPs have been associated with embryonic patterning [24], inhibition of meduloblastoma cell proliferation [25] and inhibition of glioma cell motility [26]. Inhibition of the Wnt pathway by Wnt5-a has also been shown to supress tumor activity in thyroid carcinoma [27].
This study has several limitations. First, the expression profiles in tissue collected and analyzed postnatally may not reflect embryonic expression. Consequently, whether the differences we observed are causes or consequences of the ectopic location of the thyroid remains to be tested. Second, even though clusters of genes involved in histone and chromatin function have repressed expression in ectopic thyroids, we have not formally excluded a role of differential histone methylation or acetylation on differential gene expression in ectopic thyroids. Third, the arrays used for the CNVs and methylome analyses have their own limits in definition and genome coverage. Lastly, the sample number is small but our preliminary findings justify testing a larger number of samples.
This study identifies interesting candidate pathways that may play important roles in the migration of the embryonic thyroid and provides a prototype approach for the study of congenital disorders difficult to explain by classical genetics.

Ethic Statement
This study was approved by the Ethics Committee of the CHU Sainte-Justine (ERB number 94). All the parents gave written informed consent.

Patients and Tissue Collection
We obtained flash-frozen samples of ectopic thyroid tissue removed from 3 girls aged 8, 10 and 15 yr, because it caused local symptoms (i.e. dysphagia). For controls, we used (i) thyroid tissue from 2 girls (aged 15 and 16 yr) and 1 boy (aged 4 yr) who were operated for a single hyperfunctioning thyroid nodule that had arisen in an orthotopic thyroid and (ii) commercially available RNA from normal thyroid when appropriate (Table S5).

Functional clusters and pathways analysis
We submitted the 1011 differentially expressed genes into gene ontology (GO) groups using the DAVID database (http://david. abcc.ncifcrf.gov) for cluster analysis according to Gene Ontology (GO) terms with medium or high classification stringency. To provide a refined analysis, the 100 validated gene were analyzed through Ingenuity Pathways Analysis (IPA; http:www.ingenuity. com), a software that also considers the level of gene expression. With either DAVID or IPA, the proportion of each gene in the submitted list is compared with the one in the whole genome to compute the P value of the Fisher's test, the enrichment scores (i.e. geometric mean of the inverse log of each P value) and the Benjamini-Hochberg multiple correction P value.

Expression Arrays
After surgical resection, the samples were immediately frozen in liquid nitrogen and stored at 270 Celsius until use. Total RNA was extracted as per manufacturer recommendations using the QIAzol kit (QIAGEN Inc., Ontario, Canada). RNA was DNasetreated to minimize DNA contamination. RNA quantity was measured by ND-1000 (Nanodrop, Wilmington, DE, USA). RNA quality was assessed by electropherograms on the Agilent 2100 Bioanalyzer. Microarray hybridization was performed on three different ectopic thyroids (two in duplicate, one in quadruplicate) and compared to RNA of thyroid tissue from a Caucasian female, age 68 y with gall bladder cancer (Ambion, #AM6872). Doublestranded cDNA was synthesized from 1 mg of total RNA, followed by production of antisense RNA containing the modified nucleotide 5-(3-aminoallyl)-UTP using the Amino Allyl Messa-geAmp TM II aRNA Amplification kit (Ambion, Texas, USA). After labeling with Cy3 or Cy5 (GE Healthcare Bio-Sciences, New Jersey, USA), sample pairs were hybridized onto Human Exonic Evidence Based Oligonucleotide HEEBO slides (Stanford Functional Genomics Facility, CA, USA). The oligonucleotide set consists of 44544 70-mer probes that were designed using a transcriptome-based annotation of exonic structure for genomic loci. Hybridizations were replicated with dye swap. Slides were scanned using a Molecular Devices 4000B Laser scanner and expression levels were quantified using GenePix Pro 6.1 image analysis software (Axon Instruments, CA, USA). Image acquisitions were performed with automatic photomultiplier gains (PMT) adjustment. Artefact-associated spots were eliminated by both visual and software-guided flags, as were spots with a signal/ background fluorescence ratio less than 2. The fluorescence values were imported into Acuity 4.0 software package (Molecular Devices, Union City, CA, USA). A non-linear locally weighted scatter plot (Lowess) normalization method applied to each individual block (print-tip option) was carried out using Acuity 4.0 software package (Molecular Devices, Union City, CA, USA) [28]. The identification of genes with significant differences in expression levels was performed using the significance analysis of microarray method (SAM one class) [29]. SAM estimates the percentage of genes identified by chance, the false discovery rate (FDR). We assessed the statistical significance of the differential expression of genes by computing a q-value (i.e. minimum FDR) for each gene (Table 1). Genes were considered to be differentially expressed when the absolute normalized fold change between ectopic thyroids and control was determined to be greater than 2.0 or less than 0.5 in at least one pair of the hybridized arrays. Full access to the primary array data is available on the GEO web site (http://www.ncbi.nlm.nih.gov/projects/geo/) under accession number GSE16804.

Quantitative Real-Time RT-PCR
Validation of the expression levels of 100 genes of interest was carried out using TaqMan low density array (TLDA) technology (Applied Biosystems, Ontario, Canada). Probes and primers have been selected with the publically available software http://www5. appliedbiosystems.com/tools/ and can be retrieved by using the assay ID reported in Tables 2, 3 and S4. The expression levels were normalized to the expression level of the 18S rRNA. Induced (n = 49) and repressed (n = 51) genes were selected from the 1011 differentially regulated genes found in the HEEBO expression microarray analysis. Total RNAs from thyroids were first treated with the DNA-free kit to remove residual contamination of genomic DNA (Ambion Inc.). DNA-free total RNA (175 ng) was subjected to reverse transcription using High-Capacity cDNA Reverse Transcription Kits (Applied Biosystems). An aliquot of the cDNA was mixed to TaqManH Gene Expression Master Mix, loaded on the TLDA plates then centrifuged for distribution of the material in the 384 wells. Gene target amplifications were performed in triplicate using the 7900HT Real-Time PCR System (Applied Biosystems). Sequence Detection System software version 2.2.2 (Applied Biosystems) was used for comparative gene expression analysis using the Ä Ct method. In a first analysis, expression levels found in the normal orthotopic thyroid from Ambion were compared to the mRNA levels present in the ectopic thyroids. Expression levels in the hot nodules were then compared to the levels found in the ectopic thyroids. For analysis, the cut-off log 2 value was 0.5. Then, to compare the results of the quantitative real-time RT-PCR with those of expression arrays, Pearson correlation was calculated with the free statistical software R [30].

Methylation Profiling by Methylated DNA Immunoprecipitation (MeDIP) and MeDIP-chip
The MeDIP-chip was performed using pairs of enriched methylated fraction (IP) and normal fraction (IN) of genomic DNA from our three ectopic thyroids and three controls (i.e. hyperfunctioning thyroid nodules). The methylated fraction of genomic DNA was enriched using the methylated DNA immunoprecipitation (MeDIP) assay [31] and interrogated on human Promoter plus CpG Island Tiling Arrays with a ChIP design for CpG islands and promoter regions (n = 28,226) from HG18 using 385,020 Probes selected from CGH probe bank with a median spacing of 101bp (Roche NimbleGen, Madison, WI).
Briefly, 4 mg MseI digested genomic DNA was immunoprecipitated with monoclonal mouse anti 5-methylcytidine antibody (New England Biolab, Pickering, Ontario and Abcam Inc.Cambridge, MA 02139). After washes and purification steps, immunoprecipitated material and a sample of input DNA were amplified using GenomePlex Complete Whole Genome Amplification (WGA) kit (Sigma-Aldrich, Saint Louis, Missouri 63103 USA). The resulting products (4 mg) were labeled, cohybridized and scanned by the NimbleGen Customer Service (Roche NimbleGen, Madison, WI). For each sample, NimbleScan detects peaks by searching for at least two probes above a P-value cutoff (2log 10 ) of 2 and peaks within 500bp are merged (gff files on GEO web site). Then, peak data were analyzed to compare the methylation profile between ectopic thyroids and orthotopic thyroids using the Loess normalized log 2 (ChIP/input) ratios with the one-way ANOVA tool of the Partek Genomics Suite (PGS) software. A FDR less than 0.1 was considered as significant. Full access to the primary array data is available on the GEO web site (http://www.ncbi.nlm.nih. gov/projects/geo/) under accession number GSE17581.

Array Comparative Genomic Hybridization and validation with quantitative real-time PCR
Array comparative genomic hybridization (aCGH) was performed using pairs of thyroid tissues and leucocytes from the three ectopic thyroids. We used the Nimblegen X1 HG18 whole genome CGH design (version 2). The 385,815 probes are distributed across the genome with a median spacing of 7073 bp. Probes are 60-mers, with a Tm target of 80 degrees. Labeling, hybridization, washing and scanning was performed by the NimbleGen Customer Service (Roche NimbleGen, Madison, WI). After normalization, the log 2 (test/reference) signals were analyzed using a circular binary segmentation algorithm (segMNT) with the PGS Software to identify somatically acquired segmental copy number changes. To call a copy number change, segMNT required a segment to span a minimum of 5 consecutive probes with a p-value threshold of 0.001 and a signal to noise ratio of 0.3. Then, reported regions were set at segMNT log 2 ratio below 0.3 or above 0.2 in all three samples with a p-value threshold of 0.01. Full access to the primary array data is available on the GEO web site (http://www.ncbi.nlm.nih.gov/projects/ geo/) under accession number GSE17463. Validation of the aCGH with quantitative real-time PCR was performed with the TaqMan technology. Identified CNVs were validated using TaqMan Gene Copy Number Assays from ABI. Probes and primers have been selected with the public available software (http://www5.appliedbiosystems.com/tools/cnv/) and can be retrieved by using the assay ID reported in Table S5. Figure S1 Reliability of the HEEBO expression array was confirmed by calculating the Pearson correlation coefficient (r = 0,86; P,2.2 e-16, n = 100 genes; ectopic thyroid (n = 3) vs normal thyroid -Ambion, #AM6872) between microarray and qRT-PCR results. Results are expresed in log2 ratio. Found at: doi:10.1371/journal.pone.0013420.s001 (0.01 MB PDF)