Novel NR2F1 variants likely disrupt DNA binding: molecular modeling in two cases, review of published cases, genotype-phenotype correlation, and phenotypic expansion of the Bosch-Boonstra-Schaaf optic atrophy syndrome.

Bosch-Boonstra-Schaaf optic atrophy syndrome (BBSOAS) is a recently described autosomal dominant disorder caused by mutations in the NR2F1 gene. There are presently 28 cases of BBSOAS described in the literature. Its common features include developmental delay, intellectual disability, hypotonia, optic nerve atrophy, attention deficit disorder, autism spectrum disorder, seizures, hearing defects, spasticity, and thinning of the corpus callosum. Here we report two unrelated probands with novel, de novo, missense variants in NR2F1 The first is a 14-yr-old male patient with hypotonia, intellectual disability, optic nerve hypoplasia, delayed bone age, short stature, and altered neurotransmitter levels on cerebrospinal fluid testing. The second is a 5-yr-old female with severe developmental delay, motor and speech delay, and repetitive motion behavior. Whole-exome sequencing identified a novel missense NR2F1 variant in each case, Cys86Phe in the DNA-binding domain in Case 1, and a Leu372Pro in the ligand-binding domain in Case 2. The presence of clinical findings compatible with BBSOAS along with structural analysis at atomic resolution using homology-based molecular modeling and molecular dynamic simulations, support the pathogenicity of these variants for BBSOAS. Short stature, abnormal CNS neurotransmitters, and macrocephaly have not been previously reported for this syndrome and may represent a phenotypic expansion of BBSOAS. A review of published cases along with new evidence from this report support genotype-phenotype correlations for this disorder.

INTRODUCTION

Bosch–Boonstra–Schaaf optic atrophy syndrome (BBSOAS) is a recently described autosomal dominant disorder characterized by optic atrophy and/or hypoplasia, developmental delay, and intellectual disability (Brown et al. 2009; Al-Kateb et al. 2013; Bosch et al. 2014; Chen et al. 2016). This disorder is caused by alterations in NR2F1, a highly conserved nuclear receptor protein that regulates transcription (MIM: 132890). NR2F1 belongs to family of orphan nuclear receptors. It has two functional domains, the DNA-binding domain (DBD) and the ligand-binding domain (LBD), both of which are highly conserved across the members of the nuclear receptor family; the ligand for the LBD is, however, still unknown. It is very similar to its mouse counterpart, COUP-TF1, which has been extensively studied using knockout models. Many physiological roles of COUP-TF1 have been delineated as a result—neurogenesis (Zhou et al. 1999; Yamaguchi et al. 2004), eye development (Tang et al. 2010), cortical patterning (Faedo et al. 2008), to name a few. Most variants in NR2F1 described thus far have been missense variants that lead to haploinsufficiency or dominant negative effects and are predominantly located in the two functional domains. A smaller number of indels and larger deletions have also been reported (Brown et al. 2009; Al-Kateb et al. 2013; Chen et al. 2016). Here we describe two unrelated patients, with novel, de novo, missense variants in NR2F1, with new clinical features, not yet described for this syndrome. We also review the phenotypes of published cases and discuss their relation to the location of the variants. Using molecular modeling we are also able to demonstrate the putative effect of the two missenses on protein function, generating a hitherto undescribed molecular model for the LBD of NR2F1 in the process.

RESULTS

Clinical Presentation

Case 1

The patient was a third child born after normal pregnancy to a 32-yr-old mother. Apgar scores were 9 and 9 at 1 and 5 min, respectively, and birth weight was 6 pounds and 14 ounces (2785 g). The patient has had problems with feeding since birth and his psychomotor development was delayed. At age 1 yr the patient was functioning at a 2-month developmental level. One episode of infantile spasms was reported at 6 months that resolved with a single dose of ACTH; prophylaxis was not begun. He underwent bilateral inguinal herniorraphy at 10 months and subsequently received a gastrostomy tube to facilitate feeding. He was evaluated in the genetics clinic at 14 yr of age. His weight was 27.60 kg, height was 147 cm, and head circumference was 55 cm. Facial features were notable for mild dolichocephaly, mild midface hypoplasia, deep-set eyes with short palpebral fissures, large protruding ears with simplified helix and deficient lobuli, bilateral skin tags on the posterior aspect of the ear lobes, and micro- and retrognathia with large appearing teeth (Fig. 1A). The patient was nonverbal with very limited social interaction and needed assistance with ambulation. Self-stimulating behavior was observed in the form of constant rocking back and forth and repetitive movements of his arms. Mild facial weakness, global hypotonia, diminished muscle mass, and motor apraxia were evident. Deep tendon reflexes were bilaterally brisk. The patient was able to localize touch and there was no cerebellar involvement. Short stature, both height and weight consistently below the third percentile, was recorded between 14 and 19 yr (Fig. 1B,C). Bone age was ∼12.5 years at a chronological age of 15 yr, representing a delay of greater than 3 SDs. Brain magnetic resonance imaging (MRI) was unremarkable and limited ophthalmologic evaluation disclosed mild optic nerve hypoplasia bilaterally. An awake and asleep computer-assisted prolonged video EEG recorded left occipital onset seizure with secondary generalization, dysrhythmia grade 3, generalized and multifocal spikes and sharp waves, and right temporal intermittent rhythmic delta activity. Cerebrospinal fluid (CSF) examination was significant for a very low level of hydroxyindoleacetic acid (5-HIAA) of 19 nmol/l (normal range 67–140 nmol/l) and a low level of homovanillic acid (HVA) of 68 nmol/l (normal level of 145–326 nmol/l), whereas concentrations of 3-methyldopa, tetrahydrobiopterin, neopterin, and 5-methyltetrahydrofolate reductase were within the laboratory reference range. The overall phenotype was significant for optic nerve hypoplasia, severe intellectual disability with absent speech and limited ambulation, hypotonia, dysmorphic facies, short stature with delayed bone age, and low levels of 5-HIAA and HVA in CSF.

Figure 1.

Phenotype of Case 1 with Bosch–Boonstra–Schaaf optic atrophy syndrome. (A) Facial features: deep-set eyes, midface hypoplasia, simplified ear lobules, protruding ears, micrognathia, and retrognathia. (B) Height chart and (C) weight chart. Yellow crosses indicate serial measurements of the patient's height and weight. All marks fall consistently below the third percentile.

Case 2

The patient was the first child born at 35 weeks via C-section following premature rupture of membranes to a 25-yr-old mother. The pregnancy was complicated by hyperemesis gravidarum and maternal hypertension in the last trimester. Birth weight was 5 pounds and 15 ounces (2336 g). Feeding difficulties were noted at birth and a nasogastric tube was placed for the first 3 days, after which the infant was bottle-fed. At 5 1/2 weeks a heart murmur was noted and subsequent echocardiogram revealed a bicuspid aortic valve as well as a mildly dilated aortic root. These have remained stable in follow-up echocardiograms. Hypotonia was noted around 15 months. Anisometric amblyopia in the left eye was diagnosed at 17 months, and the patient underwent left lateral rectus recession for left dissociated horizontal deviation at 5 years. There was significant developmental delay, both motor and speech, with the patient functioning at a 6–8-month-old level at 15 months and at 2-year-old level at 5 years of age. The patient exhibited significant repetitive behavior such as rocking or banging her head against a wall, particularly while trying to fall asleep. On examination at age 2 years height was at 26th centile and weight was at 64th centile (Fig. 2). The patient was mildly dysmorphic with slight hypertelorism, prominent synorphrys, simplified cupped ear helices, and large head (head circumference at 2 yr was 53 cm, 97th percentile). Wide proximal phalanges as well as mild hypotonia were also noted. The patient had normal hearing. MRI was not performed.

Figure 2.

Phenotype of Case 2 with Bosch–Boonstra–Schaaf optic atrophy syndrome. (A) Facial features: hypertelorism, prominent synorphrys, simplified cupped ear helices, and large head. (B) Height chart and (C) weight chart. Yellow crosses indicate serial measurements of the patient's height and weight.

Genomic Analyses

Genetic testing for fragile X, Prader–Willi, and Angelman syndromes (including UBE3 sequencing), MECP2, TGFβR1, and TGFβR2 were normal. PTEN gene analysis revealed a variant in the 5′ UTR c.488_486dupGGC, inherited from the father. This variant was not interpreted as being phenotypically significant, given that the father was normal. Karyotype performed on skin fibroblasts was normal (46,XY). Chromosome microarray analysis revealed two copy-number changes. The first was a maternally inherited 58-kb duplication, 3p13(71,152,133–71,210,298) ×3, which included the 5′ UTR, exon1, and exon 2 of FOXP1 gene (OMIM: 605515). Although deletions in this gene have been associated with intellectual disability, language impairment, and autism (Hamdan et al. 2010), this variant was classified as a familial variant with no phenotypic significance. The second, a paternally inherited 700-kb interstitial deletion at 7p14.3 (33,772,968–34,472,618) ×1, covers a region that includes two genes (BMPER and NPSR1). This deletion was reported as likely benign. Clinical whole-exome sequencing performed on the proband and his parents identified a novel, de novo, heterozygous, missense variant in the NR2F1 gene (NM_005654.5), c.257G>T, p.Cys86Phe (Table 1). The C86 position of NR2F1 is a highly conserved in orthologous protein sequences from human to Drosophila. No variants resulting in C86F in NR2F1 are reported in any of the major population databases (ESP [Exome Variant Server], ExAC [Lek et al. 2016], gnomAD [Lek et al. 2016], 1000 Genomes [1000 Genomes Project Consortium et al. 2015], dbSNP [Sherry et al. 2001]). A synonymous change encoding this residue is reported in ExAC at a very low frequency (3 in 45,280). Bioinformatics algorithms were applied to predict the effect of C86F (SIFT score, 0 [Kumar et al. 2009]; PolyPhen-2 score, 1.000 [Adzhubei et al. 2010]), which agreed that the variant might be damaging to the protein.

Table 1.

Novel NR2F1 variants

Gene	Transcript	Chromosome position (GRCh37/hg19)	HGVS	Type/effect	Genotype	ClinVarID	dbSNP	Allele frequency ExAC/gNOMAD (all)	In silico tools (PolyPhen-2, SIFT)	Predicted effects with molecular modeling and dynamics simulations
NR2F1	NM_005654.5	5:92920986	c.257G>T, p.Cys86Phe	Substitution/missense	Heterozygous	SCV000599463	No record	No record	Damaging	Decreases stability of the zinc-binding domain
		5:92929391	c.1115T>C, p.Leu372Pro	Substitution/missense	Heterozygous	SCV000599464	No record	No record	Damaging	Disrupts dimerization interface of the ligand-binding domain

Genetic testing for fragile X syndrome, Cornelia de Lange syndrome, and chromosomal microarray were normal. Clinical whole-exome sequencing performed on the proband and her parents found two paternally inherited missense variants, NOTCH1 (NM_01617.4) c.7115GArg2372Gln and SGSH (NM_000199.3) c.220C>T, p.Arg74Cys. A cardiac ultrasound on the father could not be performed to help differentiate if the bicuspid aortic valve is a coincidental finding in the patient or if the NOTCH1 variant is likely responsible for this cardiac finding. A novel de novo, heterozygous, missense variant in the NR2F1 gene (NM_005654.5), c.1115T>C, p.Leu372Pro was identified (Table 1). The L372 residue of NR2F1 is conserved in orthologous protein sequences from human to Caenorhabditis elegans. No variants resulting in L372P in NR2F1 are reported in any of the major population databases (ESP, ExAC, gnomAD, 1000 Genomes, dbSNP). Bioinformatics algorithms were applied to predict the effect of L372P (SIFT score, 0; PolyPhen-2 score, 1.000), which agreed that the variant might be damaging to the protein.

Molecular Modeling and Dynamics Simulations

Molecular modeling was used to provide a more detailed model to inform our clinical interpretation of these novel variants.

NR2F1 C86F Alters Zinc Ion Coordination and Decreases Stability of the Zinc-Binding Domain

The Cys86Phe substitution lies in a highly conserved zinc-finger (ZF) nuclear hormone receptor-type domain, which is the first of two ZF domains within the DNA-binding domain. Each zinc-finger domain has four cysteine residues that coordinate the bound ion, and Cys86 is among them. Previous studies have established that changes in the DNA-binding domain disrupt transcriptional activity (Chen et al. 2016), and we hypothesized that changes to Cys86 may also have this effect. Thus, a homology-based protein model was developed and used in an in silico test for the effect of C86F on the structure and dynamics of the ZF domain (Fig. 3A–C; see Methods for details). In our duplicate molecular dynamics (MD) simulations of the DNA-bound form, we identified differences in the time-dependent dynamics of the protein wherein the variant caused the structure to adopt a different conformation (Fig. 3D). We quantified this relationship using principal component (PC) analysis. Further, the distances between residues that make up the zinc-binding site are more variable for C86F compared with wild type (wt) (Fig. 3E,F). Changes in the geometry of the zinc coordination center for C86F simulations were evident (Fig. 3C), indicating that the variant is likely to alter zinc ion coordination and thereby the stability of the domain. We performed MD of apo (no DNA bound) NR2F1 in triplicate and at two temperatures (300°K, room temperature, and 360°K) in order to better understand the potential effects of C86F on stability and folding. Using the same analyses as we used to quantify differences in the DNA-bound form, significant and consistent distortions to the zinc-binding site were observed, as measured by the separation of coordinating cysteine residues (Fig. 3G). Altered dynamics were enhanced at higher temperature (Supplemental Fig. S1). At a global level, C86F induced larger root-mean-square deviation (RMSD), a measure of the difference in the global shape of NR2F1, at both temperatures, as well as increased root-mean-square fluctuation (RMSF), a measure of each residue's mobility (Supplemental Fig. S2). To understand the potential functional significance of the observed changes to NR2F1 in our simulations, we also simulated four previously established pathogenic variants: R112K, R142L, C128R, and F110del (Fig. 4). We have summarized the impact of variants using a simple metric—the separation between the center of mass (COM) of the DNA and protein. This simple metric summarizes across three types of effects: alteration in specific DNA interactions by Arg residues (R112K and R142L), local loss of structure due to Zn+2 binding site alteration (C86F and C128R), and collapse of the structure due to loss of F110 in the hydrophobic core. Although the detailed changes in the structure of NR2F1 and its interactions with DNA differ by each variant, the impact on how the protein interacts with DNA is consistent across pathogenic variants and the VUS observed in our patient case.

Figure 3.

C86F leads to distortion of the zinc-binding site. The zinc-finger domain of NR2F1 and bound DNA are shown in cartoon representation. Solvent atoms are hidden for clarity, and zinc ions are represented by purple spheres. (A) Molecular model of NR2F1 zinc-finger (ZF) domains shows the structural role of C86. (B) A representative from among the largest deviations in wt simulations is shown. Although the zinc ion has moved, the geometry of the binding site is preserved. (C) Throughout simulations of F86, the geometry of the zinc-binding site is significantly altered and zinc ion interacts with residues outside of the binding site. (D) Plotting each frame from the simulation as a point in the dominant principal component (PC) subspace, the two proteins adopted different conformations. (E,F) The distance (Å) between the Cα atom of residue 86 and of two of the other zinc-coordinating residues. Both measures showed stability for the duplicate wt simulations, but instability for C86F. (G) The distance (Å) between zinc-coordinating residues across apo simulations. Probability density plots and structural representatives (inset) show that across triplicate simulations of apo-NR2F1 at 300 K and 360 K, C86F led to greater instability around the ion-binding site.

Figure 4.

Pathogenic variants in the DBD lead to greater separation from DNA. Four previously reported pathogenic variants were simulated similarly to C86F and wt. We measured the distance between the center of mass (COM) of NR2F1 and of the bound DNA fragment. The separation between the COM of each molecule was monitored, and all pathogenic variants lead to a significant increase. However, the two variants that directly altered zinc-binding residues, C128R and C86F, exhibited the greatest effect.

NR2F1 L372P Disrupts the Dimerization Interface of the Ligand-Binding Domain

The Leu372Pro substitution lies in a ligand-binding domain (LBD) that is highly conserved (97% identical) compared with the paralog NR2F2, but for which the structure has not been experimentally determined. LBD of NR2F2 is also the dimerization domain. NR2F1 is known to bind DNA as dimers (Schrader et al. 1996; Zhang and Dufau 2001), and therefore it is very likely that the NR2F1 LBD is also the dimerization domain. We generated a molecular model of the NR2F1 LBD dimer from the experimental structure of an NR2F2 homodimer. In this model, the majority of the dimer interface is composed of hydrophobic residues within an α-helical secondary structure that we will refer to as the binding helix. L372P lies in the middle of the binding helix. Proline placed in the middle of an α-helix in soluble proteins is well known to destabilize the helix (Li et al. 1996). Consequently, we hypothesize that the observed substitution will destabilize the ligand-binding helix. To test this hypothesis, we employed MD simulation and measured the conformational stability of the binding helix. Compared with wt and across triplicate simulations, heterodimers and homodimers of L372P exhibited lower α-helix content and greater conformational variability, as measured by RMSD (Fig. 5). Therefore, MD simulations support our hypothesis that the binding helix is destabilized by L372P. Because it is known that NR2F1 binds DNA as a dimer, we hypothesize that loss of dimerization due to Leu372Pro would impair DNA binding and thereby alter transcriptional regulation.

Figure 5.

L372P destabilizes the ligand-binding domain (LBD) dimer interface. (A) Examination of our model of the LBD dimer revealed that L372 is in the center of the α-helix that makes up much of the dimer interface. (B) The dimerization helices from each monomer are arranged in parallel, leading to the L372P and G368D residues from each to be facing one another. (C) We observed loss of α-helical content within the dimerization helices after MD simulation. (D) L372P also leads to alteration of binding helix conformation as quantified by larger RMSDs.

Our structural model also provided a mechanistic hypothesis for the effect of another observed NR2F1 LBD variant identified in association with BBSOAS and located within the binding helix—Gly368Asp (Chen et al. 2016). The binding helices from each monomer are aligned in parallel. In the dimerized form, each Glycine at position 368 faces one another. Substitution to a negatively charged Asp will affect dimerization by adding a larger side chain and a negatively charged surface. These two features are likely to destabilize NR2F1 dimerization.

DISCUSSION

NR2F1 (nuclear receptor group 2, family 1), also known as COUP-TF1 in mice, belongs to the group of orphan nuclear receptors. NR2F1 has been shown to be one of the transcription factors expressed in bones of the skull and limbs where it interacts with bone morphogenetic protein-4 (BMP-4) (Feng et al. 1995) and can be responsible for aberrant mineralization levels at these sites, define susceptibility to pathological osteopenic changes, and cause fate specification. Alterations in these activities can be hypothesized to lead to the delayed bone development that we report for the first time as being associated with this syndrome, as evidenced by short stature and delay in bone maturation in Case 1. NR2F1 influences the balance of cortical patterning between motor and sensory areas in the brain, a phenomenon known as arealization (Armentano et al. 2007). Expression of NR2F1 in iPS cells has recently been shown to lead to their differentiation into cells that have the same neurotransmitter profiles as GABAergic neurons (Teratani-Ota et al. 2016). The neurotransmitter abnormalities observed in Case 1 can be hypothesized to be the result of disturbances in arealization as well as differentiation of neuronal cells. Further characterization of NR2F1 and its role in the neuronal circuitry is needed to gain more insight into these mechanisms. Functional studies have confirmed a dominant negative effect (Chen et al. 2016) for previously reported pathogenic missense variants in the DNA-binding domain (DBD) of NR2F1. Case 2 exhibits a slightly minor phenotype of the disease, as do previously reported LBD variants. The previously reported cases (Chen et al. 2016) had only mild ophthalmologic findings, as does our case. No hypotonia was noted in the previous cases, whereas in our case hypotonia was noted in infancy but has resolved with growth. An MRI was not performed in our patient; however, the previously reported cases have normal MRIs. Including the two patients in this report, a total of 30 patients with BBSOAS have been described in the literature (Table 2). Of these, there are eight patients with whole-gene deletions of NR2F1 in Chr 5 encompassing NR2F1 among other genes, ranging in size from 582 kb to 5 Mb. Smaller indels were reported in 5 of 30 cases. Missense mutations are the most common mechanism of variation, with 17 of 30 patients having de novo mutations. Of the missense mutations, 10 of 17 lie in the DNA-binding domain, four are in the start codon, and only three are in the ligand-binding domain, making a case for the DNA-binding domain being a mutational hotspot. Although phenotypic characteristics range from mild to severe for patients with all variants, patients with variants in the ligand-binding domains seem to have a milder degree of developmental delay, no hypotonia, no speech defects, no seizures, or no repetitive behaviors. These milder features do corroborate a possible genotype–phenotype relationship for this syndrome, with variants in the ligand-binding domain leading to a milder phenotype (Table 2).

Table 2.

Variants in NR2F1 and phenotypic features of patients with Bosch–Boonstra–Schaaf optic atrophy syndrome

GRCh37(Chr 5) Sl. no.	NR2F1 Chr. position	NM_005654.5 cDNA nomenclature	AA change	Domain	Inheritance	Gender	Age	Development	CVI	Hypotonia	Oromotor dys	MRI
1	Chr 5: g.92920731	c.2T>G	p?	Affects start codon, the third codon is also ATG, so alternate start possible	De novo	Male	7 yr	+	Optic nerve hypoplasia	+	+	Normal
2	Chr 5: g.92920731	c.2T>G	p?	Affects start codon, the third codon is also ATG, so alternate start possible	De novo	Female	4 yr	+	ONA	−	+	Thin corpus callosum
3	Chr 5: g.92920731	c.2T>C	p?	Affects start codon, the third codon is also ATG, so alternate start possible	De novo	Female	3 yr	+	CVI	+	−	NA
4	Chr 5: g.92920732	c.2T>C	p?	Affects start codon, the third codon is also ATG, so alternate start possible	De novo	Female	12 yr	+	ONA	+	+	Thin corpus callosum, asymmetry in the fiber directionality in the right temporal, temporo-occipital, and parietal lobes when compared with the contralateral hemisphere
5	Chr.5: g.92920731_ 92920733del	c.2_4delTGGinsG	p?		De novo	Male	43 yr	+	ONA	+	+	Bilateral occipital stenogyria, bilateral hippocampal malrotation, and abnormal configuration of both occipital horns
6	Chr.5: g.92920832_ 92920842delinsCGCC	c.103_113delinsC	Gly35Argfs*8	DNA-binding domain	De novo	Male	30 yr	+	ONA	+	+	Normal
7	Chr 5:g.92920986	c.257G>T	Cys86Phe	DNA-binding domain, zinc ion bind	De novo	Male	14 yr	+	mild b/l	+	+	Normal
8	Chr 5: g.92921020del	c.291delC	Tyr98Thrfs*2	DNA-binding domain	De novo	Male	21 yr	+	ONA	+	+	NA
9	Chr.5: g.92921057_ 92921059del	c.328_330delTTC	Phe110del	DNA-binding domain	De novo	Female	6 yr	+	ONA	+	−	Thin corpus callosum
10	Chr.5: g.92921057_ 92921059del	c.328_330delTTC	Phe110del	DNA-binding domain	De novo	Male	7 yr	+	−	+	−	Mild asymmetry of lateral ventricles
11	Chr 5: g.92921064	c.335G>A	Arg112Lys	DNA-binding domain, putative DNA-binding site	De novo	Female	35 yr	IQ52	−	NA	NA	Normal
12	Chr 5: g.92921068	c.339C>A	Ser113Arg	DNA-binding domain	De novo	Female	2 yr 4 mo	Delayed	CVI	NA	NA	Small optic disc and chiasm
13	Chr 5: g.92921073	c.344G>C	Arg115Pro	DNA-binding domain	De novo	Male	12 yr	IQ48	CVI	NA	NA	Normal
14	Chr.5: g.92921111	c.382T>C	Cys128Arg	DNA-binding domain	De novo	Male	6 yr	+	ONA	+	+	Thin corpus callosum
15	Chr.5: g.92921132	c.403C>A	Arg135Ser	DNA-binding domain, putative DNA-binding site	De novo	Female	4 yr	+	CVI	+	+	Thin corpus callosum
16	Chr.5: g.92921142	C.413G>A	Cys138Tyr	DNA-binding domain	De novo	Male	15 yr	+	ONA	−	−	Delayed myelination, tonsillar ectopia, b/l symmetric loss in the periatrial white matter with secondary dysmorphic enlargement of the atria
17	Chr.5: g.92921154	c.425G>T	Arg142Leu	DNA-binding domain, putative DNA-binding site	De novo	Female	4 yr	+	ONA	+	+	Thin corpus callosum
18	Chr.5: g.92921165	c.436T>C	Cys146Arg	DNA-binding domain	De novo	Male	11 yr	+	ONA	+	+	Thin corpus callosum
19	Chr.5: g.92921192	c.463G>A	Arg155Thr	DNA-binding domain	De novo	Male	10 yr	−	−	+	−	Normal
20	Chr 5: g.92923914	c.755T>C	Leu252Pro	Ligand-binding domain, putative c	De novo	Female	18 yr	IQ 55–65	CVI	NA	NA	Not done
21	Chr.5: g.92929379	c.1103G>A	Gly368Asp	Ligand-binding domain	De novo	Male	21 yr	+	−	−	−	Normal
22	Chr 5:g.92929391	c.1115T>C	Leu372Pro	Ligand-binding domain	De novo	Female	15 mo	+	−	+, resolved with age	+	Not done
23	Chr.5: 92,742,875-93,324,530 (hg18; NCBI Build 36.1, March 2006)	582 kb		Genes included: FLJ42709, NR2F1, FAM172A, POU5F2, and one microRNA gene, MIR2277	De novo	Male	8 yr, 3 mo	+	ONA, small optic chiasm	+	+	Normal
24	Chr 5: 92,856,299–93,054,63	0.2 Mb		Genes included: NR2F1, FAM172A (partial)	De novo	Female	8 yr	+	ONA	+	−	NA
25	Chr 5: 92845157–93679748	0.83 Mb		Genes included: NR2F1, FAM172A, KIAA0825 (partial)	Paternal	Female	24 yr	IQ 61–74	CVI			Not done
26	Chr 5: 92,910,393–93,806,93	0.9 Mb* 1st pt		Genes included: NR2F1, FAM172A, KIAA0825 (partial)	NA	Male	35 yr	+	ONA	−	−	NA
27	Chr 5: 92,910,393–93,806,93	0.9 Mb* 2nd pt		Genes included: NR2F1, FAM172A, KIAA0825(partial)	Paternal	Male	2 yr	+	ONA	+	+	Thin corpus callosum
28	Chr 5: 92,878,375–94,046,21	1.2 Mb		Genes included: NR2F1, FAM172A, KIAA0825, ANKRD31	NA	Male	37 yr	IQ 70–96	ONA	−	−	NA
29	Chr 5: 91,064,110–93,896,378	2.85 Mb		Genes included: NR2F1, FAM172A, KIAA0825, ANKRD31	De novo	Female	4 yr	Delayed	CVI	NA	NA	Normal
30	Chr 5: 90,566,268–95,580,99	5.0 Mb		Genes included: NR2F1, FAM172A, KIAA0825, ANKRD32, MCTP30	NA	Female	6 yr	+	−	+	+	Normal

AA, amino acid; ADHD, attention-deficit/hyperactivity disorder; ASD, autism spectrum disorder; b/l, bilateral; CSF, cerebrospinal fluid; CVI, cortical visual impairment; mo, months; MRI, magnetic resonance imaging; NA, not available; OCD, obsessive-compulsive disorder; ONA, optic nerve atrophy; PDD-NOS, pervasive development disorder not otherwise specified; SI. no., serial number; yr, years.

Molecular modeling predicted the effects of these variants on protein function and exemplified the practical application of these methods as recently advocated (Oliver et al. 2016). The application of molecular modeling and computational methods was used to characterize the two reported novel variants in NR2F1, the first in the DBD and the second in its LBD. Alterations to zinc-binding residues, including our novel C86F, were associated with the more significant alterations in our simulations, compared to other pathogenic variants. This may be due to two effects: First, the local unfolding of the zinc-binding site is a more drastic structural change, compared to alteration of a surface property, as in the substitution of arginine residues. Second, we have not considered the difference in contact between motif-containing DNA and nonspecific DNA. Therefore, clearer differences from wt for previously identified pathogenic variants may be achieved by considering a more specific metric, but the effect of C86F on NR2F1 structure is clear.

This report reiterates a possible genotype–phenotype correlation, with variants in the ligand-binding domain leading to milder phenotypes. It also potentially expands the phenotypic spectrum of BBSOAS to include short stature with delayed bone age and low levels of serotonin and dopamine metabolites (5-HIAA; HVA) and macrocephaly. We cautiously interpret these new features as phenotypic expansions of this newly described and currently evolving syndrome. The complete description of the phenotype will be better established as more cases of this rare disease are reported. Interestingly, although one of our patients had amblyopia, we did not see optic atrophy in either of our patients. Only 15 of the 28 previously reported cases have documented optic nerve atrophy (ONA; Table 2), whereas seven of the remaining 13 had a diagnosis of cerebral visual impairment (CVI; Table 2). Detailed ophthalmologic exams were available in 27 of 30 cases. Small or pale optic discs, strabismus, and esotropia and amblyopia were commonly encountered (Table 2). These findings suggest that that optic atrophy is an inconsistent feature. Renaming the syndrome as BBSS may be more reflective of this phenotypically evolving recently discovered syndrome associated with NR2F1 mutations.

METHODS

Sample Collection and Whole-Exome Sequencing

Whole-exome sequencing (WES) was performed on genomic DNA extracted from all samples submitted. The exome was captured utilizing a custom reagent developed by the Mayo Clinic and Agilent Technologies, targeting 19,456 genes and 187,715 exons using 637,923 probes to capture a 54.1-Mbp total region. Sequencing was performed on an Illumina HiSeq 2500 Next-Generation sequencing instrument, using HapMap Sample NA12878 as an internal control. Paired-end 101-bp reads were aligned to a modified human reference genome (GRCh37/hg19) using Novoalign (Novocraft Technologies). Sequencing quality was evaluated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). All germline variants were jointly called through GATK Haplotype Caller and GenotypeGVCF (McKenna et al. 2010). Each variant was annotated using the BioR Toolkit (Kocher et al. 2014) and subsequently evaluated for clinical relevance. Sequence coverage is detailed in Table 3.

Table 3.

Sequencing coverage

	Total mapped reads (millions)	Reads on target (%)	Duplication reads (%)	Target reads >20% (%)	c.257G>T (reads—total:variant)
Case 1
Proband	115.7	73.0	8.2	>97	197:97
Mother	130.6	72.0	9.4	>97	205:0
Father	90.1	73.0	7.4	>97	148:0
Case 2
Proband	135.6	87.4	9.3	>97	237:107
Mother	161.0	87.3	7.6	>97	319:0
Father	150	86.6	11.9	>97	262:0

Molecular Modeling and Molecular Dynamic Simulations

We generated a molecular model of the DNA-bound zinc-finger domain of NR2F1 using the homologous (54% identical) crystal structure of RXRA (1BY4 [Zhao et al. 2000]) and validated it against the experimentally solved solution structure (DOI: 10.2210/pdb2ebl/pdb.). Explicit solvent molecular dynamics (MD) simulations were carried out using NAMD (Phillips et al. 2005) and the CHARMM27 with CMAP (Mackerell et al. 2004) force field. The wt system was solvated (Grubmüller and Groll), ions added to 150 mM NaCl, and computational mutagenesis used to make the C86F model. The wt and C86F proteins were minimized, heated, and equilibrated over a combined 5 nsec each, and in duplicate (NPT). A further 25 nsec of simulation trajectory was generated for each (NVT) and the final 20 nsec analyzed. Additional and independent duplicate 2 nsec explicit solvent MD simulations were generated using the CHARMM force field within Discovery Studio (BIOVIA). Simulations of wt and C86F were also performed without DNA present (apo) and at two temperatures: 300°K and 360°K. Triplicate simulations of each condition were performed in NAMD with systems prepared similarly to above. For each, 22 nsec of simulation was generated and the final 15 nsec analyzed. Thus, a total of 384 nsec (264 nsec apo, 120 nsec holo) of MD simulation was leveraged to characterize the effects of C86F.

We additionally studied four previously reported pathogenic variants: R112K, R142L, F110del, and C128R. Triplicate explicit solvent MD simulations of these four variants, wt, and C86F were generated by a similar procedure to the above but with the following modifications: First, the bound DNA fragment was restrained in simulation to limit its dynamics. Second, the protein was initially constrained using harmonic restraints and these restraints were slowly released over 5 nsec. Finally, 20 nsec of production simulation was generated and analyzed.

Analysis was carried out using custom scripts, leveraging VMD (Humphrey et al. 1996) and the Bio3D R package (Grant et al. 2006). Protein structure visualization was performed in PyMol (“The PyMOL Molecular Graphics System. Version 1.5.0.3”) and VMD. Prior to analysis, all trajectories were aligned to the initial wt conformation using Cα atoms. PC analysis was performed using Cα atoms in Cartesian space.

ADDITIONAL INFORMATION

Data Deposition and Access

The variants mentioned in this report were deposited in ClinVar (http://www.ncbi.nlm.nih.gov/clinvar/) on 9/12/2017. The ClinVar IDs are SCV000599463 and SCV000599464 for the NM_005654.5:c.257G>T and NM_005654.5:c.1115T>C variants, respectively.

Ethics Statement

The proband and/or parents were consented for sample collection and subsequent analysis under a protocol approved by the institutional review board of the Mayo Clinic (IRB ID 12-009346: Whole-Exome Sequencing of Diagnostic Odyssey cases evaluated at the Individualized Medicine Clinic). Written informed consent was obtained from the probands’ parents for publication and accompanying images.

Acknowledgments

We thank the patients and their family for their participation in this report.

Author Contributions

C.K., M.T.Z., E.W.K., and D.B.-V. designed the study. C.K., M.T.Z., Z.N., M.J.F., R.A.U., E.W.K., and D.B.-V. gathered the data. C.K., M.T.Z., and E.W.K. analyzed the data. C.K., M.T.Z., Z.N., M.J.F., R.A.U, E.W.K., and D.B.-V. wrote the paper.

Funding

We thank the Mayo Clinic Center for Individualized Medicine (CIM) for supporting this research through the CIM Investigative and Functional Genomics Program.

Competing Interest Statement

The authors have declared no competing interest.