Somatic mutations are not observed by exome sequencing of lymphocyte DNA from monozygotic twins discordant for congenital hypothyroidism due to thyroid dysgenesis.

BACKGROUND/AIMS: Congenital primary hypothyroidism (CH) is a rare pediatric disorder estimated to occur in about 1:2,500 live births. Approximately half of these cases entail ectopic thyroid tissue, which is believed to result from a migration defect during embryogenesis. Approximately 3% of CH cases are explained by mutation(s) in known genes, most of which are transcription factors implicated in the embryology of the thyroid gland. Surprisingly, monozygotic (MZ) twins are usually discordant for CH due to thyroid dysgenesis, suggesting that most cases are not caused by transmitted genetic variation. One possible explanation is somatic mutation in genes involved in thyroid migration occurring after zygotic twinning. Such mutations should be observed only in the affected twin.
METHODS: To test the hypothesis of somatic mutation, we performed whole exome sequencing of DNA from three pairs of MZ twins discordant for CH with ectopic glands.
RESULTS: We found no somatic mutations exclusive to any of the three affected twins or in any of the unaffected twins.
CONCLUSION: Either somatic mutations are not significant for the etiology of CH or else such mutations lie outside regions of the genome accessible by exome sequencing technology.

Introduction

Congenital primary hypothyroidism (CH) is a relatively rare disorder estimated to occur in about 1: 2,500 live births [1] . This pathology is classified principally into two groups: (1) dysgenesis (85%) defined by an embryologic defect in tissue development and (2) dyshormonogenesis (15%) with anatomically normal tissue development. CH due to thyroid dysgenesis (CHTD) occurs with an incidence of 1: 4,000 live births, either with an observed ectopic thyroid gland (75% of cases) or with complete athyreosis (25% of the cases) [2] . Thyroid gland ectopy results from a defect in the migration of precursor cells from their origin at the back of the tongue to their final location at the base of neck, while athyreosis (total absence of the thyroid) results from a lack of differentiation of the thyroid during the embryogenesis. It is unclear whether these two groups are mechanistically related or result from different embryogenic mechanisms. CHTD is predominantly nonsyndromic and sporadic (98% of cases are nonfamilial) [3], and has a female and ethnic predominance [4]. Mutation(s) in transcription factors implicated in the morphogenesis of the thyroid gland (NKX2.5, NKX2.1, PAX8 and FOXE1) have been found in only 3% of CHTD cases, mostly sporadic ones [1, 3–5]. Linkage studies have excluded these genes in the rare examples of multiplex families with CHTD [6]. All these points, together with the high discordance rate (92%) between monozygotic (MZ) twins [7] , argues against inheritance of highly penetrant germline mutations as major causal factors. In contrast, somatic mutations of high penetrance arising during embryogenesis (after the stage of MZ twinning) are potential explanatory factors.

Whole exome sequencing (WES) is a new and powerful tool for identifying genetic causes of human diseases. This technique involves sequencing the annotated protein coding exons of the genome with nearby flanking intronic regions (or splicing sites, defining by intron-exon junctions) [8]. A major proportion of known mutations causing high-penetrance genetic disorders are found in or adjacent to the coding regions, so WES is a good approach to explore potential genetic causes of CH. The observed rate of somatic mutation, as assayed by whole genome sequencing of DNA from multiple tissues from the same individual, is extremely low [9, 10]. Indeed, two recent studies have failed to detect genetic or epigenetic differences between MZ twins discordant for different pathologies [11, 12]. Nonetheless, somatic mutations remain of interest as possible causal factors for rare non-Mendelian disorders.

This study aims to identify potential genetics causes for CHTD, under the hypothesis that high penetrance somatic mutations in protein coding regions of the genome are responsible for the disorder. To test this hypothesis, we performed WES in three pairs of MZ twins discordant for CHTD.

Patients and Methods

Patients

Three MZ twin pairs (2 female pairs and 1 male pair; fig. 1) discordant for CHTD were studied. Peripheral blood was obtained by venipuncture, and DNA was prepared from peripheral leukocytes using standard methods. This study was approved by the Ethics Committee of the CHU Sainte-Justine. All the parents and participants provided written informed consent.

Fig. 1

Pedigrees of the 3 pairs of twins. Squares = Males; circles = females; filled symbols = affected individuals; open symbols = unaffected individuals.

Exome and Sanger Sequencing

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. Sequencing was performed at the genomic platform of the McGill University and Genome Quebec Innovation Centre. Primary data in the form of fastq files of raw sequence plus quality scores were processed by the bioinformatics team of the Réseau de médecine génétique appliquée (RMGA). The analysis pipeline for short read alignment to the human consensus genome assembly (v19), variant calling and functional annotation included GATK (Genome Analysis Toolkit), SAMTools and ANNOVAR, using the Integrative Genome Viewer (IGV) to visualize results. Raw sequences were independently analyzed using NextGene (SoftGenetics, Inc.).

The output from bioinformatics analysis provided data for more than 50 variables for each individual, including mutation location, read coverage of wild type and variant alleles in the exome samples, allele frequencies in public databases (1000 Genomes, NHLBI Exome Variant Server), functional predictions of pathogenicity for amino acid missense variants, and the genotype quality (GQ) score generated by the GATK Unified Genotyper algorithm [13]. We sought to define stringency levels for various quantitative metrics to optimize false-positive versus false-negative genotype calls, to control for both common variants (with high minor allele frequencies) and to exclude likely exome-generated technical artifacts. We reduced this complexity through hierarchization in order to generate a testable number (n) of variables [GQ, exonic + splicing annotation, 1000 Genome minor allele frequency, control exome minor allele frequency, variant zygosity in project exome samples, variant read coverage, variant read percentage, and gene families known to be prone to exomic artifacts (mucins, HLAs, MAGEs, NBPFs, PRAMEs), etc.]. A key component was anonymized access to results from more than 1,000 local control exomes of individuals who were either healthy or had pathologies unrelated to thyroid function, generated using the same sequencing instrumentation technology and bioinformatics analysis pipeline. We generated all possible combinations of the reduced set of these n variables taken k at a time, with each variable combination yielding a set of a candidate variants and genes. Based on the reduction of the number of variants and the overlapping results, we chose five variables that generated the best set of candidate rare variants in or near protein-coding regions: GQ score, exonic or adjacent intronic location, minor allele frequency of ≤ 1% in 1000 Genomes data, variant present in 5 or fewer local control exomes (MAF 0.5%), and variant present in only one member of each twin pair (fig. 2). Only somatic variants specific to the affected twin were of biological relevance, but we looked for variants specific to the unaffected twin in each pair as well, to help define the total amount of somatic mutation.

Fig. 2

Breakdown of the experimental protocol including the bioinformatics pipeline and filtering of WES data.

Candidate variants were retested by PCR-based Sanger sequencing. PCR amplicon primers were designed with Primer3Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi/). Sanger fluorescent sequencing with capillary electrophoresis was performed at the MUGQIC using standard procedures, and sequence chromatograms were analyzed using MutationSurveyor (Soft Genetics, Inc.).

Results

Clinical Characteristics of Patients and Twins

The clinical characteristics of the 3 pairs of twins are reported in table 1. The following findings are of note and consistent with a previous report [7]: (1) all affected twins have a relatively modest increase in TSH levels at screening when compared to the confirmatory TSH values at diagnosis, which is evidence of fetal blood mixing between the healthy (T4 ‘donor’) and affected (T4 ‘receiver’) twins, and (2) on average, the birth weight of the affected twin is similar to that of his or her healthy sibling, another argument against environmental fetal distress as a risk for CHTD [2].

Table 1

Clinical characteristics of the affected and healthy MZ twins

Twins	Sex	Birth weight, g	At neonatal screening		At scintigraphy			Diagnosis
			TSH (<15), mU/l	Total T4 (120–350), nmol/l	TSH (0.4–11), mU/l	Free T4 (11–23), pmol/l	Tg (<20), μg/l
J1	F	2,080	38	120	257.75	4.38	1.9	athyreosis1
J2	F	2,140	0	180	NA2	NA2	NA2	healthy
J3	M	2,783	22	71	>324	3.42	34.5	apparent athyreosis1
J4	M	2,500	0	143	1.35	15.8	57.90	healthy1
J7	F	2,460	28	74	338.98	2.44	NA	thyroid ectopy
J8	F	2,020	<2	NA	NA3	NA3	NA3	healthy

Values in parentheses represent the normal range. Tg = Thyroglobulin; NA = not available.

Also confirmed by thyroid echography.

Normal thyroid function at 2 months of age (TSH 2.20 mU/l and free T4 12.32 pmol/l).

Normal thyroid function at 6 months of age (TSH 2.47 mU/l and free T4 12.43 pmol/l).

WES Analysis

For WES, we typically obtain 100- to 120-fold median coverage of targeted exonic regions using the Agilent SureSelect plus Illumina HiSeq technologies. Following short read alignment, variant calling and gene annotation by our bioinformatics analysis pipeline, each exome typically yields more than 100,000 candidate variants, of which approximately 25,000 are either in or immediately adjacent to coding exons. By direct comparison of exonic variants without additional filtering, our twin exomes contained approximately 400 candidate variants specific to one or the other twin. As this is far larger than can be reasonably expected based on known somatic mutation rates [9], it was clear that a substantial proportion of these represented technical artifacts, either due to the exome sequencing itself, or introduced during the multiple cycles of PCR during exome library prep, or due to incorrect read alignments to the genome consensus (for example, in segmental duplicate regions or in genes with processed or unprocessed pseudogenes). Therefore, we included additional restrictive filters to obtain a manageable number of candidate somatic variants for verification by PCR-based Sanger sequencing. Based on a comprehensive comparison of potential filters and quantitative metrics we required a minimum GQ score (generated by the GATK Unified Genotyper algorithm) metric of 90 for all genotype calls. This reduced the number of candidate variants which are different between twins to fewer than 50. We used our set of in-house control exomes to further eliminate variants arising recurrently through technical artifacts especially misalignments (and at the same time, to eliminate common variants potentially arising through recurrent mutation at hypermutable sites, such variants being of unlikely biological relevance). Restricting to variants occurring in fewer than 5 such control exomes, there were fewer than 10 potential somatic variants in each twin. The short read genomic alignments for these were each visually inspected with IGV, to remove additional likely artifacts, such as variants arising uniquely at the 3′-most position in all variant-containing reads, and variants seen in genomic neighborhoods with multiple nearby rare variants (indicative of problematic short read alignment). After applying all of these filters, each twin (both the affected and the unaffected one) had typically fewer than 5 candidate somatic variants. These variants were retested by PCR-based Sanger sequencing. None of the tested variants in the 3 twin pairs were validated by this resequencing. All candidate variants were found to be artifacts, either false positive (no variant seen in either twin member) or false negative (variant seen in both twins) ( tables 2, 3; online suppl. fig. 1–3; for all online suppl. material, see www.karger.com/doi/10.1159/000365393). False positives of high genotype quality were presumed to be due to PCR-induced mutations during the exome library prep; false negatives were presumed to be due to stochastic sampling bias in the twin exome lacking the variant.

Table 2

Numbers of variants obtained after different filters and validation

Samples	Total detected variants	≥90	1000 Genomes frequency (≤0.01 + blank)	1000 Genomes frequency + exonic or splicing	1000 Genomes or splicing +external control ≤5	Variants potentially discordant after filtering	Discordant variants validated
J1/J2	135,161	72,335	13,769	2,975	437	8	0
J3/J4	142,487	81,545	16,411	3,353	416	2	0
J7/J8	139,213	79,953	15,807	3,135	418	4	0

Details of numbers of variants obtained after each step of filtration and the number of variants validated for the 3 twin pairs are shown.

Table 3

List of variants named after filtering data

Sample	w/variant by exome	Gene	Chrom.	Position	Nucleotide changed	Acid amine changed	Variant in IGV (Jx/Jx)
J1	+	STX5	11	62592936	NM_003164: exon6:c.C499GNM_001244666: exon6:c.C499G	NM_003164: exon6:p.H167D|NM_001244666: exon6:p.H167D	−/−
J1	+	FRG1B	20	29625985	ncRNA_splicing:FRG1B:NR_003579(NR_003579: exon4:c.528+1G>A)		−/−
J1	+	MAML3	4	140811125	NM_018717: exon2:c.C1465A	NM_018717: exon2:p.Q489K	−/−
J2	+	AOC2|AOC3	17	41003601	NM_003734: exon1:c.C241Gdownstream: AOC2:NM_009590|downstream: AOC2:NM_001158	NM_003734: exon1:p.P81Adownstream: AOC2:NM_009590downstream: AOC2:NM_001158	−/−
J2	+	ANKRD20A8P	2	95513877	ncRNA_exonic: ANKRD20A8P:NR_003366		−/−
J2	+	ANKRD20A8P	2	95513880	ncRNA_exonic: ANKRD20A8P:NR_003366		−/−
J2	+	ANKRD20A8P	2	95513889	ncRNA_exonic: ANKRD20A8P:NR_003366		−/−
J2	+	EPPK1	8	144940540	NM_031308: exon1:c.C6882T	NM_031308: exon1:p.G2294G	−/−
J3	+	GANAB	11	62400952	NM_198335: exon7:c.G661CNM_198334: exon6:c.G595C	NM_198335: exon7:p.A221PNM_198334: exon6:p.A199P	−/−
J4	+	JAKMIP1	4	6086690	NM_144720: exon5:c.835_836delNM_001099433: exon5:c.835_836del	NM_144720: exon5:p.279_279delNM_001099433:exon5:p.279_279del	−/−
J8	+	LOC441666	10	42832211	ncRNA_exonic: LOC441666:NR_024380		−/−
J8	+	WASH2P	2	114355998	ncRNA_exonic: WASH2P:NR_024077		−/−
J8	+	SDHA	5	256514	NM_004168: exon15:c.G1974C	NM_004168: exon15.p.P658P	−/−
J8	+	TCERG1	5	145838680	NM_006706: exon4:c.A672GNM_001040006: exon4:c.A672G	NM_006706: exon4:p.Q224QNM_001040006: exon4:p.Q224Q	−/−

All variants (name of gene, number of chromosome, position of the variant, change of nucleotide and acid amine) for each twin pair are shown. It is indicated in which twin the variant is viewed with IGV.

Discussion

To our knowledge, this study is the first analysis to search for genetic variants in MZ twins discordant for CHTD. However, our study shows no genetic difference between such discordant MZ twins. Our results are not unexpected given that no genetic differences were reported between MZ twins in other similar studies, either healthy or discordant for different medical conditions (i.e., multiple sclerosis, VACTERL association-type congenital malformations, and renal agenesis) [10, 11, 14].

At least three possible explanations may account for our finding. First, WES covers only 1–2% of the genome and does not assess variation in noncoding gene regulatory elements except when these are immediately adjacent to coding exons or in most of 5′ - or 3′-nontranslated mRNA sequences. Moreover, a small fraction of bona fide protein-coding exons are refractory to either hybrid capture or next-generation sequencing with current technologies. Somatic mutations discordant in our twins and accounting for CHTD might be found in regions outside the accessible coding sequences of the genome. Second, epigenetic differences have been shown between MZ twins [15], which could account for CHTD [16]. Indeed, genetically identical twins show differences in their patterns of DNA methylation and an increased acetylation of histones during their lives [15, 16]. However, our previous studies did not reveal any epigenetic differences in CpG islands of promoter regions between normal eutopic and dysgenetic ectopic thyroid tissues [17, 18]. The next logical step is therefore to analyze the epigenetic profile outside promoter regions and CpG islands.

Finally, MZ twins might share the same variants in the heterozygous state, the difference lying in the possible random monoallelic expression of these variants in thyroid tissues of the affected twin [19]. Monoallelic expression of TPO gene has been reported in CH due to dyshormonogenesis [20]. This hypothesis is difficult to test because it requires normal and ectopic tissues from MZ twin pairs to look for biased allelic expression, whereas obtaining thyroid biopsies from healthy individuals is problematic.

In conclusion, we have shown that MZ twins discordant for CHTD have an identical protein-coding genome as assessed by WES. This implies that molecular causes of CHTD might be either due to genomic changes outside the accessible regions of the genome or due to yet unidentified epigenetic differences. Consistent with our findings, it has recently been shown that the rate of mutation across the whole genome is very low between healthy MZ twins [21]. It raises also the possibility that biased somatic monoallelic expression of heterozygous variants might explain phenotypic discordance in otherwise genetically identical MZ twins. This latter hypothesis would also be consistent with the occurrence of familial cases in 2% of CHTD [3] and would imply that the rate of random monoallelic expression in thyroid tissue is also about 3–5%. Given that the CHTD incidence is about 1: 4,000 live births [2], it implies that rare heterozygous variants in gene involved in thyroid development occurs in 1: 1,000 live births. Then, if we postulate that each gene in our genome has equal chances to be mutated, it would suggest that a pool of at least 20 genes are susceptibility markers for CHTD. This new stochastic hypothesis of the cause of CHTD is consistent with the polygenic hypothesis derived from animal models [22] and is also in keeping with the mainly sporadic occurrence of CHTD and the discordance rate of 92% observed in MZ twins. Finally, given that Caucasians have concentrated lethal variants in their genome [23], the stochastic hypothesis might also account for the enrichment of CHTD cases in Caucasian populations when compared to African populations [4]; however, this hypothesis, as many others, does not provide any obvious explanation for the female predominance observed in CHTD. Nevertheless and even with its limitations, this stochastic monoallelic expression hypothesis provides a new model in which the load of heterozygous variants within a functional genomic module should be considered as potentially disease-causing [24] as long as proper functional studies validate the impact of these variants on thyroid morphogenesis.