Demonstration of Autosomal Monoallelic Expression in Thyroid Tissue Assessed by Whole-Exome and Bulk RNA Sequencing

Abstract

Background:

Congenital hypothyroidism due to thyroid dysgenesis (CHTD) is a disorder with a prevalence of 1/4000 live births, the cause of which remains unknown. The most common diagnostic category is thyroid ectopy, which occurs in up to 80% of CHTD cases. CHTD is predominantly not inherited and has a high discordance rate (>92%) between monozygotic (MZ) twins. The sporadic nature of CHTD might be explained by somatic events such as autosomal monoallelic expression (AME), given that genes expressed in a monoallelic way are more vulnerable to otherwise benign monoallelict genetic or epigenetic mutations.

Objective:

The aim of this study was to search for complete (90%) AME in normal and dysgenetic thyroid tissues.

Methods:

Aggregated analysis of whole-exome and bulk RNA sequencing was performed on two ectopic thyroids, four normal thyroids, and the human thyroid cell line Nthy-ori.

Results:

A median of 5062 (range 2081–5270) genes per sample showed sufficient numbers of heterozygous single nucleotide polymorphisms to be informative. The median monoallelic expression represented 22 (range 16–32) of the informative genes for each thyroid sample. Examples of genes displaying AME are FCGBP, ZNF331, USP10, BCLAF1, and some HLA genes; these genes are involved in epithelial–mesenchymal transition, cell migration, cancer, and immunity.

Conclusions:

AME may account for the high discordance rate observed between MZ twins and for the sporadic nature of CHTD. These findings also have implications for other pathologies, including cancers and autoimmune disorders of the thyroid.

Introduction

Congenital primary hypothyroidism (CH) is a relatively rare disorder estimated to occur in around 1/2500 live births. This condition is classified principally in two categories: first, dysgenesis defined as a defect during embryologic development; and second, dyshormonogenesis, a functional disorder within an anatomically normally located thyroid. Permanent congenital hypothyroidism due to thyroid dysgenesis (CHTD) occurs with an incidence of 1/4000 live births, usually with an ectopic thyroid gland (75%) or athyreosis (15%) (1). The ectopic thyroid results from an arrest in the migration of the precursor cells to their final location at the base of neck, while athyreosis (total absence of thyroid) results from a lack of differentiation or disappearance of the thyroid during embryogenesis. At most, 3% of CHTD cases are explained by a mutation in transcription factors implicated in the embryology of the thyroid gland (NKX2.5, NKX2.1, PAX8, FOXE1, and GLIS3) (2). Indeed, CHTD is predominantly non-syndromic and sporadic (i.e., 98% of cases are non-familial) (3), has a female and ethnic predominance (4), and has a discordance rate of 92% between monozygotic (MZ) twins (5), which implies that somatic epigenetic and/or genetic mutational events might account for this discordance. However, significant somatic methylation profile differences were not observed between ectopic and orthotopic thyroids (6), nor were somatic mutations found by exome sequencing of lymphocytic DNA from MZ twins discordant for CHTD (7). This said, other epigenetic causes might account for differences in the phenotypic expression of the disease. Widespread autosomal monoallelic expression (AME) with initial non-predetermined (“random”) choice between maternal and paternal alleles defines an unusual class of disease-causing genes (8,9). Recently, biased allele-specific expression has been observed in at least 3% of the genes expressed in thyroid tissue (10), but the exact proportion of genes subject to a complete monoallelic expression in human thyroid remains unclear. Genes expressed in a monoallelic way show decreased expression (11) and may lead otherwise benign monoallelic mutations to cause loss of function in a proportion of cells (12). Consequently, differences in autosomal somatic gene expression might account for thyroid diseases including dysgenesis (13), cancer, and autoimmune disorders. In this context, this study assessed whether thyroid tissue (ectopic and orthotopic) displays genes subject to AME (Fig. 1).

FIG. 1.

Schematic representation of the autosomal monoallelic expression assessed by whole-exome and bulk RNA sequencing. Color images available online at www.liebertpub.com/thy

Patients and Methods

Patients and tissue samples

Flash-frozen samples of ectopic thyroid tissue were obtained that had been removed from two girls (aged 10 and 8 years) due to local symptoms (dysphagia). These two girls were treated with levothyroxine and had normal thyroid hormones levels before surgery (Supplementary Table S1; Supplementary Data are available online at www.liebertpub.com/thy). For controls, healthy thyroid tissue from four children was used (three girls aged 8, 12, and 16 years, and one boy aged five years), who were either prophylactically operated (two RET-mutated individuals of known MEN2 families) or operated for papillary thyroid cancer (Supplementary Table S1). The human Nthy-ori cell line was also examined (14). Finally, to assess for somatic mutations in ectopic thyroids, whole-exome sequencing (WES) of leucocyte DNA of cases #1 and #2 was performed and compared with matched data from thyroid samples of these two cases. This study was approved by the institutional research ethics committee, and parents gave informed consent.

DNA and RNA extraction from tissue samples

DNA was extracted using a pureLink genomic DNA Mini kit (Life Technologies, Burlington, Canada). RNA was extracted using a miniRNeasy kit (Qiagen, Toronto, Canada).

WES and bulk RNA sequencing

Sequencing was performed in the McGill University and Genome Quebec Innovation Centre. WES was performed using the Agilent 50 Mb SureSelect exon capture library, followed by 2 × 100 nt paired-end sequencing on the Illumina HiSeq 2000 instrument as described previously (7,15). Paired-end and strand-specific RNA-sequencing was performed in five controls and two case samples (see sample characterization in Supplementary Table S1).

Analysis of bulk RNA sequencing

Analyses were carried out with aligned RNA-sequencing data, where alignment was based on the established TopHat 2 pipeline. To assess allelic expression, high-quality variant calls of exome sequencing data (generated separately), the variant call files were used for local genotype imputation/phasing (in RefSeq genes) using IMPUTE2 software with 1000 Genomes reference haplotypes (Phase I, vs3) in two main steps: pre-phasing and imputation. Variants with an r ² value of <0.8 and phase genotypes were filtered out using SHAPEIT2 with 1000 Genomes reference haplotypes. In the case of singleton single nucleotide polymorphisms (SNPs), an attempt was made to phase based on common heterozygous sites in paired-end reads. Allele-specific analyses for transcriptome utilized a careful realignment process to phased genomes, coupled with combined tests for allelic abundance of each transcript as recently described (van de Geijn, BioRxiv, 2014; https://dx-doi-org.web.bisu.edu.cn/10.1101/011221). To assess the thyroid sampling quality, cluster dendrograms were established by unsupervised hierarchical clustering of RNA sequencing data using the software Clustal W (www.clustal.org/clustal2/; Supplementary Fig. S1).

Combined analysis of WES and bulk RNA sequencing

To define monoallelic expression, variants heterozygous in WES and nominally or seemingly homozygous in RNA sequencing data were identified. To generate a list of testable candidate variants, the number of genes was reduced using a hierarchization pipeline for each sample (Fig. 2). First, to select genotypes from the exome data, the software SAMtools was used (16), and the variants were filtered through the following criteria: artefacts outside the exome were filtered out, number of reads in exome data (≥10), and genotype quality score (QUAL <50). Only good-quality variant calls (genotypes) were used in the latter analysis. Second, to assess which genes showed complete AME, RNA-sequencing data of the informative genes were filtered through the following criteria: percentage of reads in RNA (either <10% or >90% of the altered or reference reads, respectively), inclusion of genes with two or more variants (to eliminate false positives), and exclusion of imprinted genes (as listed in www.geneimprint.org). Candidate variants were retested by polymerase chain reaction (PCR)-based Sanger sequencing of copy and genomic DNA. PCR amplicon primers were designed with Primer3Plus (www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi/).

FIG. 2.

Breakdown of the experimental protocol, including bioinformatics pipeline and filtering of whole-exome and bulk RNA sequencing data.

Comparison of WES data from matched DNA leukocytes and ectopic thyroids

For WES, 100–120-fold median coverage of targeted exonic regions was typically obtained using the Agilent SureSelect plus Illumina HiSeq technologies. Following read alignment, variant calling, and gene annotation by the bioinformatics analysis pipeline, each exome typically yielded >100,000 candidate variants, of which approximately 25,000 are either in or immediately adjacent to coding exons. To assess for somatic mutations in ectopic thyroids, the genotype of all variants called in the RNA sequencing data (after filtrations) was verified by comparing the number (and percentage) of reads between thyroids and leukocytes.

Results

Gene expression clustering of thyroid tissues

To assess the quality of the thyroid tissue samples with respect to contamination with other tissue types, an unsupervised hierarchical clustering of the RNA sequencing data was performed. As compared to other human cells (fibroblasts, monocytes, and B- and T-lymphocytes), the gene expression profiles derived from all thyroid tissues clustered together and were distinct from other human cells in a cluster dendrogram (Supplementary Fig. S1, left side branch). These results indicate that the thyroid RNA samples were clean.

Exome sequencing, RNA sequencing, and autosomal monoallelic expression in thyroid tissue

To determine whether thyroid tissue shows autosomal monoallelic expression (AME), RNA-sequencing combined with exome sequencing was performed on seven different thyroid samples: 9078 genes with good-quality exonic SNPs were covered by WES. To exclude non-informative variants (homozygous variants) and false-positive variants, a hierarchization pipeline was used for each sample (Fig. 2). First, a median of 5062 (range 2081–5270) genes per sample showed sufficient numbers of heterozygous SNPs to be informative. Then, the median complete AME represented 22 (range 16–32) of the informative genes for each thyroid sample, and AME represented a small fraction of the informative genes (0.46%; range 0.31–0.96%). Numbers of specific (non-recurrent) genes in ectopic thyroids are comparable to that of control tissues. Next, validation with Sanger sequencing was performed on cDNA of each sample and confirmed AME for the following genes: ZNF331, CTSF, ADRB2, USP10, BCLAF1, and FCGBP (Fig. 3). Finally, 30 genes (e.g., some HLA region genes) were recurrently (in two or more samples) subject to AME, and seven of these genes are located on chromosome 6 (Table 1). Only a few genes belong to the same pathways (assessed in KEGG pathway database; www.genome.jp/kegg/pathway.html): MAP2K3 of the MAPK pathway showed AME in all samples; DTX2 (AME in control 1) and CTBP2 (AME in cases 1 and 2) belong to the Wnt and Notch pathways. Of note, given that only two ectopic thyroids are available for analysis, no formal statistics were performed to compare gene occurrence by ontology group and/or pathways.

FIG. 3.

Monoallelic expression in autosomal genes: confirmation by Sanger sequencing of monoallelic expression in (A) normal orthotopic tissues and Nthy-ori and (B) ectopic tissue. Color images available online at www.liebertpub.com/thy

Table 1.

List of Genes Subject to AME for Each Sample: Controls (Ctrl 1–5) and Cases 1 and 2

	Ctrl 1	Ctrl 2	Ctrl 3	Ctrl 4	Ctrl 5	Case 1	Case 2
30 different recurrent genes	FRG1	FRG1	FRG1	FRG1	FRG1	FRG1	FRG1
	HLA-B	HLA-B	HLA-B	HLA-B	HLA-B	HLA-B	HLA-B
	MAP2K3	MAP2K3	MAP2K3	MAP2K3	MAP2K3	MAP2K3	MAP2K3
	PDE4DIP	PDE4DIP	PDE4DIP	PDE4DIP	PDE4DIP	PDE4DIP	PDE4DIP
	BCLAF1	BCLAF1 ^*	BCLAF1 ^*	BCLAF1	BCLAF1	BCLAF1	—
	CDK11B	—	CDK11B	CDK11B	CDK11B	CDK11B	CDK11B
	PRIM2	—	PRIM2	PRIM2	PRIM2	PRIM2	PRIM2
	—	HLA-A	HLA-A	HLA-A	HLA-A	—	HLA-A
	KANSL1	KANSL1	—	KANSL1	KANSL1	KANSL1	—
	ZNF717	—	ZNF717	ZNF717	ZNF717	—	ZNF717
	CES1	—	CES1	CES1	—	—	CES1
	TIMM23	TIMM23	TIMM23	TIMM23	—	—	—
	USP10^*	—	—	USP10^*	USP10	—	USP10
	NCOR1	—	—	—	—	NCOR1	NCOR1
	PKHD1L1	—	—	PKHD1L1	—	—	PKHD1L1
	ZNF676	—	ZNF676	—	—	—	ZNF676
	—	—	—	CCHCR1	—	CCHCR1	—
	—	—	—	—	—	CDC27	CDC27
	—	—	—	—	—	CTBP2	CTBP2
	—	—	ERAP2	ERAP2	—	—	—
	GSTM5	—	—	GSTM5	—	—	—
	—	—	—	HLA-C	—	HLA-C	—
	—	—	—	HLA-DRB1	—	—	HLA-DRB1
	HSPB7	—	—	—	—	—	HSPB7
	—	—	—	—	—	MST1L	MST1L
	NEB	—	—	—	—	NEB	—
	—	—	SGK223	SGK223	—	—	—
	—	SIRPA	—	—	SIRPA	—	—
	—	—	ZNF584	—	—	—	ZNF584
	ZNF773	—	—	ZNF773	—	—	—
47 non-recurrent genes	CLCNKA	SFTPA1	ZNF880	GBP3	CTSF^*	HLA-DQB1	SEC14L3
	OR2W3	CYP4B1	BBS12	SH3YL1	ADRB2^*	PRRC2A	ZNF331^*
	BTN3A2	DTX2	TIAM2		GPRIN2	HLA-DQA1	ZDBF2
	ANKRD65	FDFT1	WBSCR27		IGFN1	ECHDC3	MUC20
	ABO	SULT1A2			ANKRD36	TCF3	PLEKHG4B
	SNX19	CCDC57			KCNJ12	ZNF443	PSMD13
		ZNF750			TTYH2	FCGBP^*	NOMO1
					NOTCH2		FAM118A
							ANKRD36B
							ICK
							KMT2C
							PABPC1
							LILRB2
# Recurrent genes (%)	19 (76%)	9 (56.3%)	15 (75%)	21 (91.3%)	12 (60%)	15 (68.2%)	19 (59.4%)
# Non-recurrent genes (%)	6 (24%)	7 (43.7%)	4 (25%)	2 (0.7%)	8 (40%)	7 (31.8%)	13 (40.6%)
Total gene number	25	16	19	23	20	22	32

Genes located on chromosome 6 are noted in bold italic. Validated genes are identified with an asterisk.

Absence of somatic mutations in the two ectopic thyroid samples

To assess whether an early post-zygotic mutation was observed in ectopic thyroids, the exome of the matched leukocytes and the thyroid of cases 1 and 2 (leukocytes and thyroid compared for each case separately) were sequenced and compared. Although these exome comparisons initially yielded a large number of candidate somatic variants, after bioinformatics quality filtering and Sanger resequencing of some candidates, no thyroid-specific mutations were observed; in particular, there were no tissue-specific variants in any of the genes currently known to be involved in thyroid development and/or associated with thyroid dysgenesis. As expected, all genes showing AME were found to have heterozygous variants in both leukocytes and thyroid tissue from the same patient.

Discussion

In humans, genes are generally biallelically expressed, with each parental homologous gene contributing to mRNA expression. However, in some cases, only one copy is expressed, which defines monoallelic expression. There are three known general mechanisms causing monoallelic expression: (i) imprinting, in which imprinted genes are expressed preferentially or exclusively from either the paternal or maternal allele (sometimes with tissue specificity); (ii) random but clonally inherited X-chromosome inactivation in female cells; and (iii) non-predetermined (“random”) AME.

In fact, the first genes discovered to be subject to AME were those encoding immunoglobulins and T-cell receptors, by yet a different mechanism involving DNA rearrangement (17,18). However, as far as we know, this mechanism is specific to these gene families. In 1994, AME was reported in olfactory receptors genes (a family of >1000 genes), which suggested that autosomal AME might affect additional autosomal genes (19). Eventually, in 2007, a genome-wide survey analyzing clonal human blood cell lines revealed widespread AME across the autosomes (8). Over the last years, two forms of autosomal AME have been defined (20). On one hand, fixed AME implies that the allele-specific expression is conserved in daughter cells after division. Thus, fixed AME is mitotically conserved and can be studied in cell lines and tissues (or subtissue sections). On the other hand, dynamic AME is observed only during a window of time, is not mitotically transmitted, and can be studied only at the single-cell level. The present study therefore assessed fixed AME in thyroid tissues and in a thyroid cell line. Of note, given that (i) some classes of genes seem to be more prone to AME (8,12,21) and (ii) herein a recurrence of AME for a median of 15 genes per thyroid sample was observed, AME probably reflects a yet unknown tissue-specific genetic and epigenetic phenomenon (22), but is unlikely to be a completely “random” event. The somatic interplay between genetic and epigenetic mechanisms that control establishment and maintenance of widespread AME remains to be discovered (9).

In thyroid tissue, monoallelic expression of TPO has been shown in cases of thyroid dyshormonogenesis (23,24), and it has been hypothesized that PAX8 monoallelic expression in the thyroid could explain a familial case of congenital hypothyroidism associated with athyreosis (25). Moreover, RNA sequencing analysis recently showed a widespread biased allele-specific expression (an attenuated form of AME) in thyroid tissue (10). This prompted the search for autosomal genes displaying AME in thyroid tissues by using paired exome and RNA sequencing of seven different thyroid tissues.

The results show AME for a median of 22 genes per sample in ectopic and orthotopic thyroid glands, as well as in the Nthy-ori cell line. These findings confirm that genes in the thyroid are subject to AME. The number of autosomal genes subject to AME per tissue was relatively modest. This can be explained by the design of this study and by some technical limitations. First, only genes with an almost exclusive monoallelic expression were selected: the expressed allele is the almost exclusively transcriptionally active allele (>90% of reads in RNA sequencing), while the other one is transcriptionally silent. Second, a significant expression level in the thyroid (a minimum of 10 reads) was required in order to call a gene subject to AME (Fig. 2). Therefore, the results may represent only the extreme end of a continuum, beginning with allele-biased expression (or allele-specific expression (10)) and ending with high-confidence genes subject to complete monoallelic expression (20). Of note, the median of 0.46% of genes subject to AME was consistent with the 0.7% of AME observed in neural progenitor cells (26). The number of informative genes was relatively low, but this was expected because most of the SNPs are present in non-coding regions of the genome, which are not covered by exome and bulk RNA sequencing. For example, the median number of 5062 informative genes per sample was consistent with the genome-wide study of AME conducted in B-cells by Gimelbrant et al., which showed about 3900 informative genes per sample (8).

CHTD is predominantly a sporadic (random) event, and ectopic thyroid, the main form of CHTD, results from an arrest in the migration of the precursor cells. By exome sequencing, it was previously shown that MZ twins discordant for thyroid dysgenesis do not differ in the DNA extracted from nucleated blood cells (7), making early post-zygotic somatic mutation a less likely cause of CHTD. Similarly, no somatic epimutations were observed in the ectopic thyroids (6). The absence of any validated thyroid-specific somatic mutations in the present study is consistent with all these previous findings. Taken together, these results imply that early post-zygotic coding region mutations, or altered epigenetic markings, are unlikely in CHTD. In contrast, the extensive AME observed in all thyroid samples assessed in the present study might well account for the discordance between MZ twins and the sporadicity observed in CHTD. Indeed, genes subject to AME might contribute to cellular (clonal) diversity (8) and are more vulnerable to otherwise benign heterozygous genetic or epigenetic mutations (27).

Only three recurrent genes are common in the two ectopic glands (CDC27, CTBP2, and MST1L), and they do not belong to similar cellular pathways (as assessed by the KEGG pathway database; www.genome.jp/kegg/pathway.html). Of note, AME genes found in other studies are enriched for transmembrane receptors (8) and glycoproteins, which are essential to cell adhesion and organ development (12). Similar examples of genes displaying AME in the ectopic thyroid tissues are ZNF331 and FCGBP; both are involved in cell division and cell migration (28,29). ZNF331 encodes a zinc finger protein containing a KRAB (Krüppel-associated box) domain found in transcriptional repressors. Recently, it has been shown that ZNF331 expression was reduced in human esophageal cancer by methylation of the promoter and that cell migration is suppressed by ZNF331 expression (28). FCGBP (the Fc fragment of the IgG binding protein) is a surface protein involved in the epithelial–mesenchymal transition and in tumor progression in gallbladder cancer (29). FCGBP expression is reduced in papillary and follicular thyroid cancer, which might suggest a role in thyroid cell division and migration (30). Remarkably, it was found that seven genes recurrently subject to AME in the thyroid were located on chromosome 6 (HLA-A, -B, and -C; DRB1; BCLAF1; PRIM2; and CCHCR1), consistent with long-range asynchronous replication and monoallelic expression under the control of a cis-acting gene named Asynchronous replication and Autosomal RNA on chromosome 6, or ASAR6 (31). Donley et al. (31) showed that the asynchronous replication of ASAR6 was coordinated in cis with two loci on the short arm of chromosome 6, the HLA locus, and an olfactory gene cluster. However, it is acknowledged that the highly polymorphic nature of MHC Class I and II genes may impact alignment efficiency of RNA-seq reads and lead to spurious allelic imbalances (32). Consequently, the frequent monoallelic events in these genes require additional validation. This suggests that even if the choice of the monoallelically expressed allele is a stochastic (non-predetermined) event, the location of the genes subject to AME is probably determined by genetic and/or epigenetic marks across the genome (9).

Stochastic variable gene expression among genetically identical cells has been implicated in processes including determination of cell fates (33,34) and the development of genetic diseases (11). Reduction of the stochastic expression noise, either by haploinsufficiency or by AME, increases the likelihood of impaired gene function and disease (11). A dynamic AME and transcription dynamic has been recently uncovered in embryonic cells (35). It is therefore likely that stochastic dynamic AME has even more fundamental implications for the penetrance and severity of congenital disorders, including CHTD (35). Somatic variation in the proportion of cells with normal biallelic expression compared to cell with AME of a pool of genes (9) might also lead to the broad range of phenotypes observed in CHTD, ranging from athyreosis to thyroid ectopy and thyroid hemiagenesis.

In conclusion, this study is the first to show clearly that complete AME is observed for groups of genes in ectopic and orthotopic thyroids. These findings, together with the reported tendency of surface proteins to be AME prone, suggest that a pathogenic variant expressed monoallelically in ectopic thyroid tissue might account for the almost universal discordance between MZ twins for this malformation. This latter point remains speculative, as this study is limited by the low number of ectopic thyroids, which are difficult to obtain because patients with thyroid ectopy almost never require surgery to remove the ectopic gland. This would imply that CHTD is a genetically driven stochastic event, in which morphogenic genes might interfere with cell motion, local physical chemoattractant gradients, and even local physical constraints, as suggested six decades ago by Turing (36,37). It is likely that the role of AME in sporadic congenital disorders is greater than currently appreciated (27).

Footnotes

Acknowledgments

We thank the patients and their parents for their cooperation. This work was supported by a grant from the Canadian Institutes of Health Research (MOP-130390 to J.D.). Research in pediatric thyroid diseases at Centre Hospitalier Universitaire (CHU) Sainte-Justine is supported by the Girafonds/Fondation du CHU Sainte-Justine (to J.D. and G.V.V.).

Author Disclosure Statement

The authors have nothing to disclose.

References

Deladoey

, Ruel

, Giguere

, Van Vliet

. 2011. Is the incidence of congenital hypothyroidism really increasing? A 20-year retrospective population-based study in Quebec. J Clin Endocrinol Metab, 96:2422–2429.

Narumi

, Muroya

, Asakura

, Adachi

, Hasegawa

. 2010. Transcription factor mutations and congenital hypothyroidism: systematic genetic screening of a population-based cohort of Japanese patients. J Clin Endocrinol Metab, 95:1981–1985.

Castanet

, Lyonnet

, Bonaiti-Pellie

, Polak

, Czernichow

, Leger

. 2000. Familial forms of thyroid dysgenesis among infants with congenital hypothyroidism. N Engl J Med, 343:441–442.

Stoppa-Vaucher

, Van Vliet

, Deladoey

. 2011. Variation by ethnicity in the prevalence of congenital hypothyroidism due to thyroid dysgenesis. Thyroid, 21:13–18.

Perry

, Heinrichs

, Bourdoux

, Khoury

, Szots

, Dussault

, Vassart

, Van Vliet

. 2002. Discordance of monozygotic twins for thyroid dysgenesis: implications for screening and for molecular pathophysiology. J Clin Endocrinol Metab, 87:4072–4077.

Abu-Khudir

, Paquette

, Lefort

, Libert

, Chanoine

, Vassart

, Deladoey

. 2010. Transcriptome, methylome and genomic variations analysis of ectopic thyroid glands. PloS One, 5:e13420.

Magne

, Serpa

, Van Vliet

, Samuels

, Deladoey

. 2015. Somatic mutations are not observed by exome sequencing of lymphocyte DNA from monozygotic twins discordant for congenital hypothyroidism due to thyroid dysgenesis. Horm Res Paediatr, 83:79–85.

Gimelbrant

, Hutchinson

, Thompson

, Chess

. 2007. Widespread monoallelic expression on human autosomes. Science, 318:1136–1140.

Savova

, Vigneau

, Gimelbrant

. 2013. Autosomal monoallelic expression: genetics of epigenetic diversity?. Curr Opin Genet Dev, 23:642–648.

10.

2015 Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science, 348:648–660.

11.

Cook

, Gerber

, Tapscott

. 1998. Modeling stochastic gene expression: implications for haploinsufficiency. Proc Natl Acad Sci U S A, 95:15641–15646.

12.

Gendrel

, Attia

, Chen

, Diabangouaya

, Servant

, Barillot

, Heard

. 2014. Developmental dynamics and disease potential of random monoallelic gene expression. Dev Cell, 28:366–380.

13.

Vassart

, Dumont

. 2005. Thyroid dysgenesis: multigenic or epigenetic … or both?. Endocrinology, 146:5035–5037.

14.

Lemoine

, Mayall

, Jones

, Sheer

, McDermid

, Kendall-Taylor

, Wynford-Thomas

. 1989. Characterisation of human thyroid epithelial cells immortalised in vitro by simian virus 40 DNA transfection. Br J Cancer, 60:897–903.

15.

Samuels

, Gallo-Payet

, Pinard

, Hasselmann

, Magne

, Patry

, Chouinard

, Schwartzentruber

, Rene

, Sawyer

, Bouvier

, Djemli

, Delvin

, Huot

, Eugene

, Deal

, Van Vliet

, Majewski

, Deladoey

. 2013. Bioinactive ACTH causing glucocorticoid deficiency. J Clin Endocrinol Metab, 98:736–742.

16.

, Handsaker

, Wysoker

, Fennell

, Ruan

, Homer

, Marth

, Abecasis

, Durbin

. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25:2078–2079.

17.

Pernis

, Chiappino

, Kelus

, Gell

. 1965. Cellular localization of immunoglobulins with different allotypic specificities in rabbit lymphoid tissues. J Exp Med, 122:853–876.

18.

Hozumi

, Tonegawa

. 1976. Evidence for somatic rearrangement of immunoglobulin genes coding for variable and constant regions. Proc Natl Acad Sci U S A, 73:3628–3632.

19.

Chess

, Simon

, Cedar

, Axel

. 1994. Allelic inactivation regulates olfactory receptor gene expression. Cell, 78:823–834.

20.

Reinius

, Sandberg

. 2015. Random monoallelic expression of autosomal genes: stochastic transcription and allele-level regulation. Nat Rev Genet, 16:653–664.

21.

Zwemer

, Zak

, Thompson

, Kirby

, Daly

, Chess

, Gimelbrant

. 2012. Autosomal monoallelic expression in the mouse. Genome Biol, 13:R10.

22.

Zhang

, Li

, Gao

, Egli

, Xie

, Deng

, Li

, Lee

, Aach

, Leproust

, Eggan

, Church

. 2009. Digital RNA allelotyping reveals tissue-specific and allele-specific gene expression in human. Nat Methods, 6:613–618.

23.

Fugazzola

, Cerutti

, Mannavola

, Vannucchi

, Fallini

, Persani

, Beck-Peccoz

. 2003. Monoallelic expression of mutant thyroid peroxidase allele causing total iodide organification defect. J Clin Endocrinol Metab, 88:3264–3271.

24.

Neves

, Mezalira

, Dias

, Chagas

, Viana

, Targovnik

, Knobel

, Medeiros-Neto

, Rubio

. 2010. Monoallelic thyroid peroxidase gene mutation in a patient with congenital hypothyroidism with total iodide organification defect. Arq Bras Endocrinol Metabol, 54:732–737.

25.

Hermanns

, Grasberger

, Refetoff

, Pohlenz

. 2011. Mutations in the NKX2.5 gene and the PAX8 promoter in a girl with thyroid dysgenesis. J Clin Endocrinol Metab, 96:E977–981.

26.

Eckersley-Maslin

, Thybert

, Bergmann

, Marioni

, Flicek

, Spector

. 2014. Random monoallelic gene expression increases upon embryonic stem cell differentiation. Dev Cell, 28:351–365.

27.

Chess

. 2012. Mechanisms and consequences of widespread random monoallelic expression. Nat Rev Genet, 13:421–428.

28.

Jiang

, Linghu

, Zhan

, Han

, Guo

. 2015. Methylation of ZNF331 promotes cell invasion and migration in human esophageal cancer. Curr Protein Pept Sci, 16:322–328.

29.

Xiong

, Wen

, Miao

, Yang

. 2014. NT5E and FcGBP as key regulators of TGF-1-induced epithelial-mesenchymal transition (EMT) are associated with tumor progression and survival of patients with gallbladder cancer. Cell Tissue Res, 355:365–374.

30.

O'Donovan

, Fischer

, Abdo

, Simon

, Peter

, Gerber

, Buergi

, Marti

. 2002. Differential expression of IgG Fc binding protein (FcgammaBP) in human normal thyroid tissue, thyroid adenomas and thyroid carcinomas. J Endocrinol, 174:517–524.

31.

Donley

, Stoffregen

, Smith

, Montagna

, Thayer

. 2013. Asynchronous replication, mono-allelic expression, and long range Cis-effects of ASAR6. PLoS Genet, 9:e1003423.

32.

van de Geijn

, McVicker

, Gilad

, Pritchard

. 2015. WASP: allele-specific software for robust molecular quantitative trait locus discovery. Nat Methods, 12:1061–1063.

33.

Blake

, Balazsi

, Kohanski

, Isaacs

, Murphy

, Kuang

, Cantor

, Walt

, Collins

. 2006. Phenotypic consequences of promoter-mediated transcriptional noise. Mol Cell, 24:853–865.

34.

Maamar

, Raj

, Dubnau

. 2007. Noise in gene expression determines cell fate in Bacillus subtilis . Science, 317:526–529.

35.

Deng

, Ramskold

, Reinius

, Sandberg

. 2014. Single-cell RNA-seq reveals dynamic, random monoallelic gene expression in mammalian cells. Science, 343:193–196.

36.

Turing

. 1990. The chemical basis of morphogenesis. 1953. Bull Math Biol, 52:153–197; discussion 119–152.

37.

Tompkins

, Li

, Girabawe

, Heymann

, Ermentrout

, Epstein

, Fraden

. 2014. Testing Turing's theory of morphogenesis in chemical cells. Proc Natl Acad Sci U S A, 111:4397–4402.