Determination of RET Sequence Variation in an MEN2 Unaffected Cohort Using Multiple-Sample Pooling and Next-Generation Sequencing

Multisample, nonindexed pooling combined with next-generation sequencing (NGS) was used to discover RET proto-oncogene sequence variation within a cohort known to be unaffected by multiple endocrine neoplasia type 2 (MEN2). DNA samples (113 Caucasians, 23 persons of other ethnicities) were amplified for RET intron 9 to intron 16 and then divided into 5 pools of <30 samples each before library prep and NGS. Two controls were included in this study, a single sample and a pool of 50 samples that had been previously sequenced by the same NGS methods. All 59 variants previously detected in the 50-pool control were present. Of the 61 variants detected in the unaffected cohort, 20 variants were novel changes. Several variants were validated by high-resolution melting analysis and Sanger sequencing, and their allelic frequencies correlated well with those determined by NGS. The results from this unaffected cohort will be added to the RET MEN2 database.


Introduction
Multiple endocrine neoplasia type 2 (MEN2) is a rare autosomal dominant inherited disorder with a high lifetime risk of medullary thyroid carcinoma (MTC) [1,2]. MEN2 consists of three syndromes: familial medullary thyroid carcinoma (FMTC), MEN2A, and MEN2B [1,3]. FMTC families have only MTC. MEN2A families have MTC, with at least one individual developing pheochromocytomas, parathyroid hyperplasia, or both. MEN2B patients have MTC (with or without pheochromocytoma) and other characteristic clinical features: mucosal ganglioneuromas, GI ganglioneuromas, eye abnormalities, and skeletal abnormalities including marfanoid body habitus [4][5][6][7]. MEN2 is caused by pathogenic mutations found exclusively within the RET proto-oncogene (REarranged during Transfection). These are gain-of-function dominant mutations which are commonly heterozygous missense mutations found at specific codons within RET exons 10, 11, and 13-16 and rarely found within exons 5 and 8 [1,[8][9][10]. The medical management for the patient and potentially their family members is based on the familial RET variation, which is usually determined by Sanger sequencing [1]. Discovery of a known MEN2 pathogenic RET mutation within a family leads to screening for MTC, pheochromocytomas, or parathyroid hyperplasia, and potentially prophylactic thyroidectomy to increase survival rate for the intractable, aggressive MTC. Approximately 75-80% of MTC patients have the sporadic form of MTC (i.e., isolated, nonfamilial MTC), not MEN2 [6]. Patients with apparent sporadic MTC are always tested for an RET germline mutation, in case they actually have MEN2 and require different medical management. Although there are many well-known pathogenic RET mutations causative of MEN2, it may be difficult to know if a rare or novel germline RET variant is a pathogenic mutation (patient has MEN2) or nonpathogenic polymorphism (patient has sporadic MTC).
Interpretation of rare and novel variants will increase in importance as more people are sequenced at the exome, whole genome, or targeted gene levels. Many new changes will be found with unknown clinical significance and their presence and allele frequency within the general population is of importance to help determine pathogenicity status of a variant. Consortiums like the 1000 genome and other large sequencing projects are making great progress in understanding population sequence variation. Yet more direct studies on single genes or gene panels can yield higher sequencing read coverage and more cost-effective sequencing over a smaller genetic area. Also, a particular chosen cohort can be sequenced for a particular locus, such as in the case of this study, where a cohort that was self-reported to have no personal or family history of MEN2 or MTC was sequenced for a section of the RET protooncogene where most pathogenic MEN2 causative mutations are located. RET sequence variation detected in this MEN2 unaffected population can then be added to the MEN2 RET database [8]. This data could be used for several reasons: (1) to help interpret the pathogenicity of clinically detected RET sequence variation; (2) as a reference for any future MEN2 case studies (variant was not found in those unaffected by MEN2 disease); (3) for improved genetic test design, to avoid or minimize designing probes or primers over known RET sequence variation.
To further reduce costs of sequencing large numbers of individuals, multiple samples can be pooled (without indexing) before next-generation sequencing (NGS). This was the focus of several studies that analyzed the ability to detect true variants within nonindexed pooled sample sets [11][12][13][14][15][16]. Thirty samples (60 alleles) were the maximum pooling number indicated by our prior studies and in other reports [12,14,17,18], for reproducible and accurate singleton allele detection within the pool (a singleton is a unique allele within the pool). A pool of this size was expected to a have a singleton allele read frequency of 1.67%, and with consideration of sequencing error rates and potential variance in NGS determined variant read frequencies, singleton variants are expected to be detected above a cutoff of >1% variant reads [17,18].
In this study, 136 individuals of an MEN2 unaffected cohort were sequenced on the illumina genome analyzer utilizing laboratory and bioinformatics protocols from our previous studies for nonindexed, multiple sample pooling. The pool size was limited to less than 30, which is the previously determined optimal pooling size for accurate singleton variant detection [12,14,17,18]. In total, 61 variants were detected within the MEN2 unaffected cohort, which included 20 novel variants.

Materials and Method
2.1. Samples. Peripheral blood samples from 136 adult volunteers (113 Caucasian and 23 non-Caucasians for ethnic diversity) were collected and deidentified using University of Utah IRB protocol no. 7740. The donors for this unaffected cohort were self-described as not having a personal or family history of neither medullary thyroid carcinoma nor multiple endocrine neoplasia type 2 (MEN2). The 51 samples used as controls were deidentified according to IRB no.7275 and were Sanger sequenced for RET exons 10, 11, and 13-16, including exon/intron boundaries. The "single-sample control" did not have RET mutations causative of MEN2, while the "50 pool" control contained many samples with known MEN2 causative RET mutations. The 50 pool control was sequenced on the illumina genome analyzer several times previously [17,18].

PCR, Library
Prep, and NGS. DNA samples were amplified from RET intron 9 to intron 16 using long-range PCR technology. Amplicons were normalized by SequalPrep (Invitrogen Corp, Carlsbad, CA), quantified using Quant-iT Picogreen dsDNA kit (Invitrogen Corp), and equimolar pooled before Illumina Library Prep, utilizing previously described protocols [18]. Between 27 and 29 Caucasian samples' amplicons were combined into four separate pools (P1, P2, P3, and P4) before Illumina Library Prep and NGS. The non-Caucasian cohort's 23 samples were sequenced in a separate pool (ethnic pool). The PCR-amplified RET positions 1-9180 are positions 43608691-43617870 in reference sequence NC 000010.10 (Table 1). Two controls were also included in this study, a single sample and also a pool of 50 samples. Each pool and each control were sequenced in a separate flow cell lane on the illumina genome analyzer, using single-end read chemistry.

Data Analysis.
Sequencing image files were processed and reads aligned to the RET reference sequence with SeqMan NGen version 2.1 software (DNAstar, Madison, WI), as described previously [17,18]. Reads used were of 67 base lengths since the 3 end read positions of longer reads can have an increase in sequencing background errors, as shown in previous studies [17,19,20]. As previously described, several base quality score screening thresholds (Qthreshold) evaluated for read coverage, errors (especially for outlier errors, which could be mistaken for false positives in a pool), variant read percentage, and base quality score statistics to determine the 30 Q-threshold should be used for analysis of all data sets, which minimized errors while maintaining adequate target read coverage (data not shown) [17,18].
Excluded from analysis was a region of repeats and homopolymers that caused misalignment errors in all data sets (designated "repeat region," amplicon positions 7686 to 7720). Changes from the reference sequence were designated variants, and the variant read percentage is the NGSdetermined allele frequency. The previously developed subtractive correction method of variant detection was applied wherein the control's variant read percentages (at every position and possible variant change) are subtracted from the pooled data's variant read percentages, to yield a pooled data set without background sequencing error [17,18].

Variant Validation.
A subset of the NGS-detected RET sequence variants were validated by either high-resolution melting analysis (HRM) and/or Sanger sequencing. The HRM analysis PCR primers for RET exons 13 and 15 were described previously [21]. RET intron 9 used HRM analysis primers (5 to 3 ): forward ACA CTG CAA TGT GCG GGT CA and reverse GTC CCC CAA CAA TGC TGC CC.  Sample DNA (∼5 to 15 ng/uL final concentration) was amplified and analyzed as described previously [21], except the LightScanner 32 instrument (Idaho Technology, Inc., Salt Lake City, UT) which was used for both PCR and HRM analysis. The LightScanner parameters included Uracil-DNA glycosylase step (50 • C for 10 min); polymerase activation (95 • C for 10 min); 40 PCR cycles (denaturation at 95 • C for 1 s, annealing at 62 • C for 1 s, extension at 72 • C for 4 s); formation of amplicon heteroduplexes (95 • C for 1 s, then cool rapidly to 40 • C for 10 s with ramp rate of 20 • C/s); high-resolution melting protocol (70 to 96 • C with ramp rate of 0.3 • C/s) [21]. In order to detect samples with a homozygous RET variant that could not be distinguished from homozygous wild-type samples during HRM analysis, the same procedure was performed, except wild-type DNA (∼5 ng/uL final concentration) which was spiked into the PCR reaction [22]. If needed, Sanger sequencing was used to confirm HRM determined variant results.

Next-Generation Sequencing (NGS) of the 50-Pool and
Single-Sample Controls. The sequence of the single-sample control used in this study exactly matched the reference sequence and therefore had no true variant changes from the reference sequence, only background sequencing error (Table 1). This sample was an ideal control for error rates since any variant reads from the RET reference at each sequence position reflects the background NGS error rates, and also illumina genome analyzer sequencing has demonstrated reproducible, nonuniform, sequence-specific background error rates, read coverage, and base quality scores between lanes and runs using the same version chemistry [11, 12, 17-20, 23, 24]. This single sample controls for the sequence-specific error rates within the pooled data sets by using the subtractive correction method, as described in our previous studies [17,18]. For subtractive correction, the single-sample control's variant reads at every possible sequence position, and change from the reference sequence is subtracted from the pool's variant read percentages. This yields an estimation of the pooled data without background sequencing error rates contributing to the variant read percentages (examples in Figure 1). The single-sample control and 50-pool data were also used for selection of the 30 Q-threshold used for quality screening of the data before analysis (data not shown) [17,18].
The 50-pool control demonstrated sensitivity to detect known variants with low read percentages for this NGS run and for using the subtractive correction method (as shown in Figure 1). The 50-pool contained 100 alleles and had an expected 1% singleton variant read frequency (singleton is unique within the pool). All 59 variants previously detected in the 50-pool were present at >0.5% variant reads and at similar percentage variant read values as determined in a previous NGS run with the same library (R 2 = 0.9991, Figure 1(a)). The 50-pool data had some potential false positives around the cutoff of 0.5% variant reads but after subtractive correction with the single-sample control data, and all true variants were readily detected from the background error (Figure 1(b)).

NGS of MEN2 Unaffected Cohort.
Based on our previous work and other studies, sample pools were restricted to 30 or less samples within each pool to result in a variant read percentage above 1%, the chosen cutoff for the most accurate singleton variant detection [12,14,17]. The 113 Caucasian samples were divided into 4 pools with less than 30 samples. Caucasian pool P1 had 27 samples, P2 had 29, P3 had 28, and P4 had 29 samples, with an expected 1.85%, 1.72%, 1.77%, and 1.72% singleton variant read frequency, respectively. All pool data sets were evaluated with and without subtractive correction of sequencing background error rates using the single-sample control (Figure 1(c) and data not shown). A total of 51 variants were detected in the Caucasian MEN2 unaffected cohort with >1% variant read values, of which 23 were not found in the non-Caucasian MEN2 unaffected cohort (ethnic pool) ( Table 1). The lowest singleton variant in the Caucasian data sets was in P2 with 1.12% variant read frequency. The 23 non-Caucasian samples were in one pool (ethnic pool), with an expected 2.17% singleton allele read frequency. A total of 38 variants were detected in this ethnically diverse MEN2 unaffected cohort with >1% variant read values, of which 10 were not found in the Caucasian data sets. The lowest singleton variant in the ethnic pool had 1.79% variant read frequency. The variant read percentages for each detected variant is shown in Table 1 per pool and also summarized for the four Caucasian pools. For comparison, the NCBI dbSNP allele frequency values for detected variants are also shown. All variants detected were intronic changes, except the expected common polymorphisms found in exons 11, 13, 14, and 15. Of the total 61 variants found in the MEN2 unaffected cohort, 20 variants were novel changes, not seen in the 50pool control or in NCBI dbSNP132.

Validation of NGS-Detected Variations.
Since the 136 unaffected cohort samples had not been sequenced previously by either Sanger or NGS methods, several variant locations within three pools (total of 79 samples) were chosen for validation. The high-resolution melting (HRM) analysis method, which is a rapid, closed-tube mutation scanning assay, was chosen to genotype each individual sample for validation of NGS variant detection and the NGS determined variant allele frequency (Table 2). High-resolution melting analysis detects sequence variation within the PCR amplicon using a saturating dsDNA dye and in many cases can uniquely identify each variant based on differential melting profiles (Figure 2(a)) [25][26][27][28]. HRM assay states 100% specificity and sensitivity for detection of heterozygous variants within small amplicons (<300 bp) [29]. HRM analysis was used to detect sequence variations within a section of RET intron 9 and exons 13 and 15 (Figure 2 and data not shown). RET exons 13 and 15 were chosen since they each contain a common polymorphism present in all pools that could be detected using the previously developed HRM frequency within that pool of 29 samples). Since some homozygous variants can have similar melting profiles as the wild-type sample [22], a technique that spikes wildtype DNA into the PCR reaction to allow distinction of homozygous variants was performed on any sample that appeared wild-type after testing in the first HRM assay. This technique identified four homozygous variants in RET exon 15 and one homozygous variant in RET intron 9 (Figure 2 with wild-type DNA spiked into the PCR reaction to help differentiate homozygous variants. One sample with a homozygous variant at position 174 with ("+spike") and without wild-type DNA spiked in is shown (blue traces) compared to a 174 heterozygous (red trace) and wild-type sample (black trace).

Discussion
This paper describes RET proto-oncogene sequence variation detected in an MEN2 unaffected cohort of 136 individuals. The previous genome analyzer sequenced 50pool library [17,18] was used to control for the detection of variations with low read frequency, and all known variants were detected >0.5% variant reads. With similar error rates between genome analyzer lanes of the same run [11,12,17,19,20,23,24], the singleton variants in a less than or equal to 30-sample pool should be accurately detected above background error using our previously determined cutoff value of >1% variant read frequency. The single sample controlled for background sequencing error rates across each RET sequence position and was used for the subtractive correction method of variant detection for pools [17,18]. The majority of the MEN2 unaffected cohort were of Caucasian ethnicity, while 23 samples were non-Caucasian (ethnic pool) and were used to identify RET variants within a more ethnically diverse sample set. The 136 samples were distributed into five nonindexed pools and were sequenced in five separate flow cell lanes. Using previously described protocols for bioinformatics, subtractive correction, and variant read cutoff value of >1% [17,18], a total of 61 RET variants were detected within the MEN2 unaffected cohort. Twenty of these variants were novel, not in NCBI dbSNP 132 (which includes 1000 Genome data) or found in our previous sample pooling studies on the RET proto oncogene 8 Journal of Thyroid Research [17,18]. Many of these novel changes were specific to either the Caucasian or ethnic samples and were of low variant read frequency, so they were likely to be singleton or doubleton variants within the pools. Several variants were verified by HRM and Sanger sequencing, including the novel change (c.1760−158C>G) with the lowest variant read percentage (1.12%).
The RET MEN2 database developed by the author so far has 147 entries, of which 74 are known pathogenic mutations and 62 are variants of uncertain significance. This database has been used as a model for predictions of phenotypic severity of variants of unknown clinical significance within the RET proto oncogene [30]. RET sequence variation data for these MEN2 unaffected cohorts will be added to the MEN2 RET database [8]. This variant data will help in medical interpretation of variations found in this RET protooncogene region using methods such as Sanger sequencing or NGS (targeted to the RET gene, the whole exome, or whole genome sequencing). Any RET sequence change detected in individuals with a family history of MEN2 symptoms or where MEN2 is suspected (patient with apparent sporadic MTC or Pheochromocytoma) can be compared to the available RET MEN2 database, and also to the benign RET sequence variation present in the large cohort of unaffected individuals that was generated in this paper. This highlights the importance of clinically relevant databases with not only known pathogenic changes, but also the inclusion of known benign changes for clinical test interpretation. The MEN2 unaffected cohort's variant results can also be used in comparison to variants detected in suspected MEN2 patients for case reports. A potential problem for genetic test design is unknown variants present in the location of the PCR primers, Sanger sequencing primers, or melting analysis probes. Results from this unaffected population will help with genetic test design, so that primers, and probes will not be designed over the known RET sequence changes. This paper presents sequence variation detection methods that could be used for other genes and analysis for specific cohorts (unaffected versus affected, by different ethnicities, or those with specific symptoms of disease). The resulting data can be added to locus-specific databases to help interpret the pathogenicity of clinically detected sequence variation. These validated methods can also apply to other pooled samples (such as genetic locations for GWAS followup or for testing-specific populations) and natural pools (such as mitochondrial heteroplasmy or mixed tumor populations).