ZBTB11 dysfunction: spectrum of brain abnormalities, biochemical signature and cellular consequences.

Bi-allelic pathogenic variants in ZBTB11 have been associated with intellectual developmental disorder, autosomal recessive 69 (MRT69; OMIM 618383). We report five patients from three families with novel, bi-allelic variants in ZBTB11. We have expanded the clinical phenotype of MRT69, documenting varied severity of atrophy affecting different brain regions and described combined malonic and methylmalonic aciduria as a biochemical manifestation. As ZBTB11 encodes for a transcriptional regulator, we performeded chromatin immunoprecipitation-sequencing targeting ZBTB11 in fibroblasts from patients and controls. Chromatin immunoprecipitation-sequencing revealed binding of wild-type ZBTB11 to promoters in 238 genes, among which genes encoding proteins involved in mitochondrial functions and RNA processing are over-represented. Mutated ZBTB11 showed reduced binding to 61 of the targeted genes, indicating that the variants act as loss of function. Most of these genes are related to mitochondrial functions. Transcriptome analysis of the patient fibroblasts revealed dysregulation of mitochondrial functions. In addition, we uncovered that reduced binding of the mutated ZBTB11 to ACSF3 leads to decreased ACSF3 transcript level, explaining combined malonic and methylmalonic aciduria. Collectively, these results expand the clinical spectrum of ZBTB11-related neurological disease and give insight into the pathophysiology in which the dysfunctional ZBTB11 affect mitochondrial functions and RNA processing contributing to the neurological and biochemical phenotypes.

Introduction

In 1995, two siblings with atypical methylmalonic aciduria and progressive encephalopathy were reported,[1] while the genetic cause of the disease remained elusive. The clinical presentation in this family (Family A) included intellectual disability, microcephaly, spasticity and dystonia. Atrophy of the cerebral white matter, cerebellum and brainstem, and basal ganglia signal changes were seen on MRI. The organic aciduria was later identified as combined malonic and methylmalonic aciduria (CMAMMA), which is genetically heterogeneous. Bi-allelic variants in ACSF3, encoding a mitochondrial malonyl-CoA and methylmalonyl-CoA synthetase, cause CMAMMA (OMIM 614265).[2,3] Sequencing of ACSF3 in Family A did not identify pathogenic variants (data not shown). In addition, a large Canadian newborn screening study identified individuals with variants in ACSF3 and CMAMMA but without neurological deficits.[4] We therefore performed whole exome sequencing (WES) of both siblings and their healthy mother. We identified compound heterozygous missense variants of ZBTB11 in trans in both siblings. ZBTB11 encodes the transcriptional regulator Zinc Finger (Znf) and Broad-complex, Tramtrack, Bric-à-brac (BTB) Domain Containing 11. Bi-allelic missense variants in ZBTB11 are associated with intellectual developmental disorder, autosomal recessive 69 (MRT69) (OMIM 618383), so far described in two families with intellectual disability and brain abnormalities,[5] and in a patient (ID: UPN-0706) with global developmental delay, microcephaly, cataract mild cerebellar atrophy and methylmalonic aciduria.[6]

In addition, we identified three individuals from two families (Families B and C) also with bi-allelic variants in ZBTB11, whose clinical manifestations overlapped those in MRT69.[5,6] In this study we compiled clinical and biochemical information of these five patients to expand the phenotype of MRT69. We performed in vitro chromatin immunoprecipitation with sequencing (ChIP-seq) and RNA-sequencing (RNA-seq) studies to explore the pathophysiology of ZBTB11 dysfunction at the cellular level.

Materials and methods

Genetic studies

In Family A, we performed WES in the two affected siblings (Patients FA.II-1 and FA.II-3) and in their healthy mother (Individual FA.I-2). We prepared the samples by extracting genomic DNA from peripheral blood using the SureSelectXT Human All Exon v5 (Agilent Technologies). Samples were sequenced on an Illumina HiSeq2000 instrument (Illumina Inc.), and WES data analysis was performed as previously described.[7] For Family B, genomic DNA was extracted using a Gentra Puregene kit (Qiagen) and then phenol chloroform purified. Libraries with ∼350 base inserts were prepared from 1 μg of DNA using the KAPA HTP Library Construction Kit protocol (Kapa Biosystems, Roche). Libraries were pooled in groups of eight using equal molar amounts of input for exome capture. The exome capture was performed according to Nimblegen’s SeqCap EZ (Roche) Library + UTR protocol that covers 96 Mb. Sequencing was performed on a HiSeq 2500 (Illumina Inc.) with version 4 chemistry to generate a minimum of 42 million 126 base paired-end reads. Reads were mapped to NCBI build 37 (hg19) using the Illumina aligner ‘ELAND’ [efficient large-scale alignment of nucleotide databases (Illumina Inc.)], then aligned with Novoalign v.3.02.07 (Novocraft Technologies, Selangor, Malaysia) and processed with Samtools.[8] These alignments, stored in BAM format, have been fed to bam2mpg (http://research.nhgri.nih.gov/software/bam2mpg/index.shtml), which calls genotypes at all covered positions using a probabilistic Bayesian algorithm (most probable genotype, MPG). These genotype calls were compared against Illumina Human 1M-Quad genotype chips, and genotypes with an MPG score of 10 or that showed >99.89% concordance with SNP ChIP data were kept. Sequence bases with a phred quality score <20 (Q20) were ignored. Only reads with mapping quality >30 were included in the analysis. Samples were sequenced to sufficient coverage such that 85% of the targeted exomes were called with high-quality variant detection. The variants were annotated using Annovar (www.openbioinformatics.org/annovar/annovar_gene.html). Variants detected in dbSNP (v.137), 1000 Genomes, NHLBI 6500ESP and HGMD were annotated. Filtering and analysis of non-sense, frameshift, splicing, missense variants were performed with VarSifter.[9] Variants with an allele frequency of >0.01 were removed. Using autosomal recessive inheritance filtering, we identified compound heterozygous variants in three genes: ZBTB11, PMS2 and ZNF92. In addition, we analysed data from genes causative for methylmalonic acidemia and identified a paternally inherited MMUT variant (NM_000255.4) c.878A>G p.Gln293Arg, which is a variant of uncertain significance. In Family C, WES was performed only on the proband (Patient FC.II-2) and the exonic regions were captured using Twist Human Core Exome (Twist Bioscience) followed by sequencing on NovaSeq6000 instrument (Illumina Inc.). WES data analysis was performed with GATK, v.3.6 adhering to the best practices,[10] and the casual variants were identified through filtration strategies using public databases in addition to an ethnic-specific control database (Iranome).[11]

Segregation of the ZBTB11 variants in Families A and C was studied by Sanger sequencing, using the following primers: forward CTTCTGTCCTTGCTGTTGTGCT and reverse AAAACGCCATGTCAGAACACAT; and forward ATAGGGCCGTTCACCAGTATGT and reverse TCTCAAACTGCTCAGGTGATCC (Family A) and forward TTCCCCCTGGTTACCAACTC and reverse AATGTGTCCAGAGCAGCTCATA (Family C). In Patient FB.II-1 variants were verified through RNA-sequencing data analysis.

RNA extraction from fibroblasts

Fibroblast cultures were obtained from skin biopsies of Patients FA.II-1, FA.II-3 and FB.II-1 and six control individuals. Cells were maintained in DMEM with 10% foetal bovine serum and 1% penicillin–streptomycin. For RNA extraction, the cells were washed with 1× PBS before treatment with trypsin and centrifuging. RNA extraction was performed using Paris™ Kit (Life Technologies) following the manufacturer’s instructions. RNA concentrations of the samples were measured on a Qubit fluorometer with RNA HS Assay Kit (Thermo Fisher Scientific). RNA integrity assessed on an Agilent 2100 Bioanalyzer using RNA 6000 Nano chips (Agilent Technologies).

RNA-sequencing and data analysis

Samples were prepared for RNA-seq with the Illumina TruSeq stranded mRNA-seq sample prep (Illumina Inc.). The libraries were indexed, pooled and sequenced on an Illumina NextSeq 500 (Illumina Inc.), with 75-bp single reads.

Reads were aligned to the Ensembl gene model GRCh37 (Homo_sapiens.GRCh37.87.gtf)[12] with HISAT2 v.2.1.0.[13] Bioinformatic analyses were performed in R, using packages DESeq2 v.1.32.0 from the Bioconductor project (http://bioconductor.org).[14] Simple independent sample comparisons between affected individuals and control groups were performed with DESeq2 with the default false discovery rate (FDR) using Benjamini–Hochberg set at P < 0.05. The differentially expressed genes (DEGs) discovered were further analysed by performing enrichment analysis using the R package ClusterProfiler v.4.0.5.[15] The code used is available in Supplementary material (RNA-seq data analysis). Gene set enrichment analysis and over-representation analysis were used for interpreting gene expression data; the method derives its power by focusing on gene sets, that is, groups of genes that share common biological functions.[16]

ChIP-sequencing and data analysis

ChIP of endogenous ZBTB11 from fibroblasts was performed on 5–8 million cells grown to log phase and cross-linked with 1% formaldehyde for 8 min at room temperature and quenched in 0.125 M glycine for 5 min. Cell lysis and chromatin shearing was done using the Chromatin shearing optimization kit Low SDS (Diagenode). Sonication of the crude nuclear fraction to achieve chromatin fragmentation of ∼550 bp was conducted for 10 rounds of 30 s on/30 s off using Bioruptor Plus (Diagenode). Antibody incubation was performed overnight at 4°C with 1 μg of a polyclonal anti-ZBTB11 antibody [#A303-240A; Bethyl Laboratories; previously used in ChIP-sequencing (seq) experiments],[17] and with equal amounts of Anti-Trimethyl-Histone H3 (Lys27) and control purified rabbit IgG (#17-622; Sigma Aldrich) as positive and negative controls, respectively. ChIP with Protein A Dynabeads (#10002D; Thermo Fisher Scientific) was performed and following reverse cross-linking, DNA was isolated and purified using QIAquick PCR purification kit and eluted into 30 µl of elution buffer. A PCR mix containing Power SYBR™ Green PCR Master Mix (#4367659; Thermo Fisher Scientific) was prepared according to the manufacturer’s specifications using positive and negative control primer sets using 10 μl DNA. Quantitative PCR (qPCR) was conducted on a 7900HT Real-Time PCR System (Applied Biosystems).

DNA concentration was measured on a Qubit fluorometer with DNA HS Assay Kit (Thermo Fisher Scientific). Integrity and size of immunoprecipitated DNA were assessed on an Agilent 2100 Bioanalyzer using DNA High-Sensitivity chips (#5067-4627; Agilent Technologies). Samples were prepared for sequencing using Illumina Nextera XT DNA seq (Illumina Inc.). The resulting libraries were indexed, pooled and sequenced initially on an Illumina MiSeq using MiSeq reagent kit v.3 (150 cycles) (Illumina Inc.). The libraries were diluted according to the information in the ChIP-seq QC report obtained from the MiSeq analyses using Model-based Analysis of ChIP-seq (MACS) v.2.1.2 software. The libraries were then sequenced on a NextSeq 500 (Illumina Inc.) and run with 75 bp single reads.

Sequenced reads (75 bp) were aligned to the human reference genome (GRCh37) using BWA v.0.7.17.[18] Duplicated reads were removed, sorted and filtered with Samtools v.1.8.[8] Significantly (P adjusted value <0.05) enriched peaks for ZBTB11 transcription factor binding in each ChIP-seq dataset were identified with MACS v.2.1.2.[19] Quality metrics for ChIP-seq data were assessed using the package ChIPQC v.1.28.0 from the Bioconductor project (http://bioconductor.org). For an additional quality control step of the ChIP-seq data, strand cross correlation metrics using phantompeakqualtools was computed.[20] Samples were grouped as control group and patient group and analysed using the R package DiffBind v.3.2.7.[21] DESeq2 v.1.32.0 was used to identifying statistically significantly differentially bound genomic sites. Gene associated region annotations were obtained with Genomic Regions Enrichment of Annotations Tool (GREAT)[22] and the R package ChIPseeker v.1.28.3.[23] Within control group analysis, all regions with overlapping peaks in at least four samples were grouped into consensus peaks and analysed. Functional enrichment analysis of consensus peaks and differentially bound peaks to identify predominant biological themes among these genes were performed using R package ClusterProfiler v.4.0.5.[15] Gene Ontology (GO) enrichment analysis of biological process terms used hyper-geometric statistical test with P-values adjusted for multiple testing by FDR control using the Benjamini–Hochberg method at a threshold of 0.05. The genes and connected biological processes identified in the consensus and differentially bound peaks were subsequently visualized using the gene concept network plot (cnet plot) and heatplot functions. This analysis was performed following the over-representation analysis seen at https://github.com/hbctraining/Intro-toChIPseq/blob/master/lessons/12_functional_analysis.md.

Quantitative PCR with reverse transcription

RNA from fibroblasts was converted to cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems).

Reverse transcription-quantitative PCR amplifications were performed using SYBR Green PCR Master Mix (Applied Biosystems) on a QuantStudio™ 12 K Flex Real-Time PCR System (Applied Biosystems). Samples were run in triplicate and amplification levels were calculated according to the 2−ΔΔCT method.[24] Dissociation curve analysis revealed a single product for each primer pair. Statistical significance was evaluated using a two-tailed Student’s t-test. The following primers were used: EEFSEC forward: GAGAACACCAAGTTCCGAGGT and reverse CTCTGGAATGCCCTGTGGAG; ACSF3 forward TGCGCCTGCCAGGTTCCGTGG and reverse CCTGCTCCTACACCATCCAC; TACO1 forward ACTGAAAATGGAGAAATCCAA and reverse TGGCACTTGTGGCTACTGTT.

Measurements of the activities of the oxidative phosphorylation enzyme complexes

We cultured human skin fibroblasts in M199 (Gibco) supplemented with 20% (vol/vol) foetal calf serum. Aliquots of 10 to 15 × 106 cells were washed with ice-cold PBS, frozen in liquid nitrogen and kept at −80°C until use. For the isolation of mitochondrial-enriched fractions, the pellets were thawed at 2–4°C and suspended in 2.9 ml of ice-cold 10 mmol/l Tris-HCl 7.6. We disrupted the cells mechanically with a 5 ml glass/t Teflon Potter-Elvehjem homogenizer in melting ice. After homogenization we added 0.6 ml ice-cold 1.5 mol/l sucrose and centrifuged the homogenate (10 min at 600g and 2°C). The 600g supernatant was centrifuged again (10 min at 14 000g and 2°C). The mitochondrial pellet was resuspended in 0.5 ml of 10 mmol/l Tris, pH 7.6, frozen in 50 µl aliquots in liquid nitrogen and kept at −80°C. The activities of the oxidative phosphorylation (OXPHOS) enzyme complexes were measured spectrophotometrically as described previously.[25-28] Protein concentration was determined according to Lowry et al.[29]

Data availability

The data will be available on request. Distribution of sensitive data may be subject to restrictions.

Results

Clinical data

Patients’ legal guardians and controls included in the studies gave written informed consent in compliance with the Helsinki Declaration and with the ethical committee of the institutions involved. For Family A, the study was approved by the Regional Committee for Medical Research Ethics—South-East Norway, REK 2010/1152a. For Family B, the study was approved by the National Human Genome Research Institute Institutional Review Board as part of NIH study 04-HG-0127 ‘Clinical and Basic Investigations of Methylmalonic Acidemia and Related Disorders’. Family C was not part of a research project and was referred to the diagnostic laboratory. The consent form was prepared according to the guidelines of Iranian Biomedical Research Ethics Committee, part of the National Institute for Medical Research Development, and the patients’ approval of additional research investigation was obtained.

We obtained informed written consent from Families A, B and C to perform genetic studies and publish clinical data and photos. An overview of clinical and genetic findings in the families is shown in Table 1.

Table 1

Clinical, genetic and biochemical findings

	Patient ID
Features	FA.II-1	FA.II-3	FB.II-1	FC.II-1	FC.II-2
Country of origin	Norway	Norway	France	Iran	Iran
Sex	Male	Female	Male	Female	Female
ZBTB11 variant; protein domain	Thr890Ala/Arg912Trp; C2H2/outside domains		Ile303Phe/Arg927Ter; BTB/C2H2	Arg52Trp; outside domains
Birth weight (g)	3060	3790	NA	2600	3100
Birth length (cm)	49	50	NA	42	NA
Birth OFC (centile)	2 cm <2.5th	10th	NA	1 cm <2.5th	10th
Age at last exam (years)	50	38	8	12.5	7
Height (centile)	150 cm (19 cm <3rd)	140 cm (16 cm < 3rd)	NA	146.5 cm (25th)	114.5 cm (10th)
OFC (centile)	50.5 (4 cm < 3rd)	48 (4 cm < 3rd)	NA	49 cm (<2.5th)	47.5 cm (<2.5th)
Bilateral cataracts	+	+	−	−	−
Craniofacial
Dysmorphisms	Long face and maxillary hyperplasia		+	Long face and maxillary hyperplasia
Neurological findings
Intellectual disability	+ Profound	+ Profound	+ Profound	+	+ Mild
Regression of development/speech	+/+	+/+	+/+	−/+ Dysarthria	+/+ Dysarthria
Ataxia	+	+	+	+	+
Drooling	+	+	NA	+	+
Muscle tone	Spasticity	Spasticity	Hypotonia axial	Hypotonia upper limbs	Hypotonia upper limbs
Dystonia	+	+	NA	NA	NA
Brain MRI findings
Cerebellar atrophy	+	+	NA	−	−
Corpus callosum/white matter atrophy	+/+	+/+	NA	+/−	+/−
Ventricles enlarged	+	+	NA	−	−
Molar tooth-like sign	−	+	NA	−	−
Biochemical findings at last examination
CMAMMA	+	+	+	−	−
Urine Cr (mmol/l) Reference interval NR[a]	26.9	25.1	NA	NR	NR
Urine MMA/Cr (µmol/mmol) Reference interval <2[30]	88.8	81.8	112.0	0.48	1.04
Urine MA/Cr (µmol/mmol) Reference interval <5[30]	12.2	11.3	47.5	0.81	0.94
Plasma MMA (µmol/l) Reference interval <0.30[b]	8.3	10.6	NA	NA	NA

C2H2 = cystein 2 histidine 2; Cr = creatinine; MA = malonic acid; MMA = methylmalonic acid; NA = not available; NR = not relevant.

Not specified because it depends on the hydration status of the individual.

Own reference interval.

Family A

Family A (Fig. 1A), previously reported by Strømme et al.,[1] includes two affected siblings, a male (Patient FA.II-1, Fig. 1B) born in 1970 and a female (Patient FA.II-3, Fig. 1C) born in 1982. Their healthy parents were non-consanguineous and of Norwegian descent. Patient FA.II-1 (Fig. 1B) had normal birth parameters, apart from a query of moderate birth asphyxia and a small occipitofrontal circumference (OFC). Early development seemed normal, as he was able to lift his head and gave eye contact at 6 weeks. At 15 months he could walk a few steps with support, but from this age neurological regression and growth retardation became apparent. Between 3 and 4 years of age, he was no longer able to crawl, chew or speak, and athetosis, spasticity and bilateral cataracts were noted. When last evaluated in 2020 (aged 50 years), he was nonambulatory and confined to a wheelchair. Purposeful limb movements were absent. Despite pronounced spasticity and rigidity, he displayed multiple motor stereotypies and some dystonic movements. He had severe contractures in his knees, hips, elbows and hands. A fixed full-length kyphosis of his spine prevented him to lie down flat (Fig. 1B). His plantar responses were extensor. Vision was uncertain, but hearing was good. He was able to swallow, but had profuse drooling. His growth was grossly impaired. Height was 150 cm (19 cm < 3rd centile at 19 years) (growth chart in use for Family A[31]), weight 37 kg (BMI 16.4); OFC was 50.5 cm (4 cm < 3rd centile) (growth chart in use for Family A[32]), compatible with profound microcephaly. His neuro-cranium seemed proportionally smaller than his viscero-cranium with receding forehead, deep-set eyes and prominent lower jaw. Despite a severe neurological phenotype, his general condition over the years had been stable with an astonishingly good mood and a remarkable absence of intercurrent febrile illnesses. He had unexplained osteoporosis leading to fracture of his left humerus after minor trauma. Dual energy X-ray absorptiometry scanning of lumbar vertebrae L2–L4 showed a T-score of −4.1 (osteoporosis is defined as a T score <−2.5). His osteoporosis was treated with zoledronic acid with modest improvement.

Figure 1

Pedigree, clinical and genetic findings in Families A and C. (A) Pedigree of Family A and segregation of the ZBTB11 variants. (B) Photos of Patient FA.II-1 at 50 years with dystonic movements and severe contractures in his elbows and hands. A fixed full-length kyphosis of his spine prevented him to lie down flat. (C) Patient FA.II-3 at 33 years displaying dystonia and choreoathetosis. She had multiple joint contractures, particularly in the hips, knees, elbows and hands, but did not have kyphosis of the spine. (D–M) Brain MRI examination shows marked generalized neurodegenerative changes in both affected siblings. Top row: Patient FA.II-1, bottom row: Patient FA.II.3. D and I are T1 midline sagittal images, E–G and J–L are axial T2 images: E and J are at the level of corona radiata, F and K at the level of the basal ganglia and G and L at the level of the pontomesencephalic junction. (H and M) Diffusion tensor imaging with colour-coded fractional anisotropy map corresponding to the level of the pontomesencephalic junction: D and I show marked thinning of corpus callosum and the brainstem and marked cerebellar atrophy with enlargement of the fourth ventricle and cisterna magna. (E and J) Loss of white matter and grossly enlarged lateral ventricles as a consequence of tissue atrophy. In E, there is generalized widening of the sulci indicating cortical atrophy. (F and K) Volume loss of the thalami and basal ganglia with widening of the sylvian fissures, particularly in F (in this image the angulation of the axial plane deviated somewhat from standard procedure due to fixed neck kyphosis). (G) Widening of the interfoliar and subarachnoidal spaces and unremarkable superior cerebellar peduncles (short arrow), whereas the superior cerebellar peduncles in L are elongated and thinned (long arrow) in a molar tooth-like fashion. (H and M) The molar tooth-like appearance is not present in H but readily appreciated in M (thin arrow indicates elongated and thinned superior cerebellar peduncles). Both images confirm atrophy of white matter tracts. Colour code: red indicates side to side directed fibres, blue indicates cranio-caudal directed fibres and green indicates anterior-posterior directed fibres. In H, there is normal decussation of the superior cerebellar peduncles indicated by the red colour (vertical arrow). In M, there was also normal decussation of the superior cerebellar peduncles (not shown). (N) Pedigree of Family C and segregation of the ZBTB11 variant. (O) Photograph of Patient FC.II-2 at the age of 7 years shows a prominent nose, broad nasal ridge, long philtrum and long face. (P and Q) Brain MRI images of Patient FC.II-2 at the age of 6 years. C axial T1-weighted image shows normal structure of the cerebellum with no vermis hypoplasia and molar tooth sign. Sagittal T2-weighted image shows the small size of the splenium of the corpus callosum. (R) Photograph of Patient FC.II-1 at the age of 12 years shows a prominent nose, broad nasal ridge, long philtrum and long face. (S) Sanger sequencing in Family A showing the missense variant c.2668A>G present in the patients and their healthy father and the missense variant c.2734C>T detected in the patients and their healthy mother. R = A or G; Y = C or T. (T) Sanger sequencing documented the missense variant c.154C>T in heterozygosity in the parents FC.I-1 and I-2 and in homozygosity in Patients FC.II-1 and II-2. Y = C or T.

Neuroimaging

Cerebral CT in 1983 (age 13 years) showed increased amounts of fluid over the convexities and widening of the Sylvian fissure, ventricular dilatation due to atrophy of central white matter, and partial agenesis of the vermis cerebelli. Brain MRI (Fig. 1D–G) in 2020 showed general loss of cortical grey matter, loss of white matter, enlargement of the ventricular system, reduced volume of the thalami, basal ganglia and brainstem, thinning of the corpus callosum, enlargement of the cisterna magna and cerebellar atrophy. Diffusion tensor imaging with colour-coded fractional anisotropy map (Fig. 1H) confirmed reduced volume of white substance and displayed normal decussation of the superior cerebellar peduncles at the level of the mesencephalon. Magnetic resonance spectroscopy (MRS) with short echo-time (30 TE) and intermediate echo-time (135 TE) and the voxels in the basal ganglia did not detect abnormal metabolites.

Biochemical studies

In 1983 (at 13 years of age) increased urinary excretion of methylmalonic acid was detected. No precise diagnosis was made and vitamin B12 injections and a protein-restrictive diet treatment failed to improve the condition. In 1993 (age 23 years), with better methodology, the amount of methylmalonic acid excreted in the urine was measured to be 200 µmol/mmol creatinine (normal <2). In 2012, moderately increased urinary excretion and elevated plasma concentration of methylmalonic acid were confirmed. Importantly, at this stage elevated urinary excretion of malonic acid was also recognized: 12.2 µmol/mmol creatinine (reference <5),[33] allowing for the designation CMAMMA. Values for haematology and vitamin B12, vitamin D, PTH, calcium and thyroid function were normal.

Patient FA.II-3 (Fig. 1C) was born after an unremarkable pregnancy and delivery, with normal birth parameters. Development was normal until 3 months when she appeared constantly dissatisfied with prolonged episodes of screaming. Muscular hypertonia and dystonic movements were observed at 6 months. OFC had fallen from 25th centile at birth to below the 2.5th centile. Bilateral cataracts were noted at 5 years. She had a single episode of seizures at 7 years. Neurological regression was slower than in the brother, and she managed to walk some steps with support and responded to verbal stimuli. At 10 years, neurological findings included athetosis and spasticity with extensor plantar reflexes. The OFC was 4.5 cm < 2.5th centile and height 7 cm < the 2.5th centile. EEG showed diffuse low voltage theta and delta and absence of alpha activity, and a marked response to photo stimulation.

At 33 years of age, she was able to walk some steps with support. Muscle tone was increased due to spasticity, rigidity and dystonia, and she displayed choreoathetosis as well as stereotypies and ritualistic behaviour. The plantar reflexes appeared to be extensor. She had multiple joint contractures, particularly in hips, knees, elbows and hands, but did not have kyphosis of the spine. Language was absent, but she could utter some sounds. Vision was uncertain, but hearing was good. She could swallow, but drooling was profuse. She had a receding forehead, deep-set eyes and a prominent lower jaw. She had a good mood and rarely contracted febrile illnesses. At the age of 38, height was only 140 cm (16 cm < 2.5th centile at 19 years), weight was 47.6 kg (BMI 24.3); OFC 48 cm (4 cm < 3rd centile at 20 years).

At 1.5 years brain MRI examination showed partial agenesis of the vermis cerebelli, a gracile brainstem, ventricular dilatation and atrophy of central white matter and low signals in the basal ganglia on T2-weighted images.[1] In 2019 (age 37 years), brain MRI (Fig. 1I–L) showed pathological changes similar to Patient FA.II-1 with markedly dilated ventricles, loss of white matter, very thin corpus callosum, gracile brainstem and very small vermis cerebelli. T2-weighted images showed atrophy of the basal ganglia and the thalami. There was also widening of the subarachnoid spaces. A progressive volume reduction of the brainstem and cerebellum compared to the previous MRI examination could be appreciated.[1] Axial view at the level of the mesencephalon disclosed thinning and elongation of the superior cerebellar peduncles compatible with a molar tooth-like sign. Diffusion tensor imaging with colour-coded fractional anisotropy map (Fig. 1M) showed thinning of the fibres of the superior cerebellar peduncles, which were laterally displaced due to ventricular dilatation. MRS did not reveal abnormal metabolic products, but a small increase of myoinositol of unknown significance was noted.

Biochemistry

In 1983 (age 1.5 years) increased urinary excretion of methylmalonic acid was detected similar to Patient FA.II-1. Like Patient FA.II-1, she was treated with vitamin B12 injections and protein restriction but without improvement. In 1993 (age 11 years), the amount of methylmalonic acid excreted in the urine was measured to be 300 µmol/mmol creatinine (normal <2). When reassessed in 2012 (age 30 years), elevated urinary excretion and plasma concentration of methylmalonic acid were documented. Notably, increased urinary excretion of malonic acid was also recognized: 11.3 µmol/mmol creatinine (reference <5,[33] compatible with CMAMMA as in Patient FA.II-1) (Table 1). Values for haematology and vitamin B12 and thyroid function were normal.

Family B

Patient FB.II-1 was reported as ‘subject 09’ by Sloan et al.[2] He initially presented with dysmorphic features, hypotonia and autistic behaviour and subsequently developed signs of regression, loss of speech and ambulation with occurrence of ataxia and peripheral neuropathy. He was diagnosed biochemically with CMAMMA on urine organic acid analysis at 26 months with 112 µmol/mmol creatinine methylmalonic acid and 48 µmol/mmol creatinine malonic acid (Table 1). Acylcarnitine profile and malonyl-CoA decarboxylase activity were normal. At last examination aged 8 years, the patient had profound psychomotor delay, axial hypotonia and loss of speech.

Family C

In Family C there were two affected sisters: the proband, Patient FC.II-2, and her elder sister Patient FC.II-1; their parents were first cousins (Fig. 1N). Patient FC.II-2 at 7 years and 1 month had OFC of 47.5 cm (1st centile) and height of 114.5 cm (9th centile). She had developmental delay, microcephaly, drooling, ataxia, muscle weakness, hypotonia of the upper limbs (particularly of the hands), dysarthria and chewing difficulties. At the age of 3 years, regression of motor functions of the facial muscles were noted, as the patient started to manifest drooling and dysarthria. She had intellectual disability (estimated as mild) and delayed walking. Dysmorphic features included long face, mildly malar flattening, hypertelorism, infra-orbital puffiness prominent nose, broad nasal ridge, long philtrum, short neck (Fig. 1O) and mild hyperlordosis. Brain MRI at the age of 6 years and 9 months showed thinning of the corpus callosum, white matter was otherwise normal, without cerebellar anomaly or molar tooth-like sign (Fig. 1P and Q). The proband had low blood concentration of TSH (0.12 mIU/l, reference for children 1–6 years: 0.85–6.50), and was started on Levothyroxine. Her affected sibling, Patient FC.II-1, examined at 12 years and 5 months had OFC of 49 cm (2 cm <2.5th centile) and height of 146.5 cm (25th centile). The clinical features of Patient FC.II-1 were similar to the proband Patient FC.II-2. Patient FC.II-1 manifested microcephaly, delayed walking, upper limb hypotonia (prominent in hands) and ataxia. She had chewing difficulty, and could not speak until the age of 4 years and received speech therapy from the age of 3. Speech articulation remained difficult. Patient FC.II-1 presented facial dysmorphism including long face, mild malar flattening, hypertelorism, infra-orbital puffiness, prominent nose, broad nasal bridge, long philtrum and thick lips (Fig. 1R). Urine analysis in the two sisters showed no evidence of organic aciduria, as malonic and methylmalonic acid concentration was within the normal range.

Genetic studies identified bi-allelic variants in ZBTB11

We performed WES of Patient FA.II-1, Patient FA.II-3 and their healthy mother Individual FA.I-2 (WES metrics available in Supplementary Table 1). Analysis of the WES data accordingto a recessive mode of inheritance identified in the two siblings compound heterozygous variants in the gene ZBTB11. ZBTB11 encodes a transcriptional regulator characterized by 12 Cys2His2 (C2H2) Znf domains, which mediate DNA binding, and by the BTB domain, which allows interactions with cell-specific proteins.[34] The variants identified in ZBTB11 (GenBank: NM_014415.3) were the following: Chr3(GRCh37):g.101370504T>C; c.2668A>G; p.Thr890Ala and Chr3 (GRCh37):g.101370438G>A; c.2734C>T; p.Arg912Trp. The p.Thr890Ala was located in Znf domain 11, while p.Arg912Trp was located between Znf 11 and 12 (Fig. 2A). The variants had a combined annotation dependent depletion (CADD) phred score of 26.2 and 34, respectively. Sanger sequencing verified the segregation of the ZBTB11 variants with the disease in the family (Fig. 1S and T). The ZBTB11 c.2668A>G variant (rs1278779161) is reported with an allele frequency of 3.981e-6 in the Genome Aggregation Database (gnomAD), with no homozygotes reported. The ZBTB11 c.2734C > T was not reported in gnomAD.

Figure 2

(A) Schematic drawing of ZBTB11 (O95625) with functional domains and position of the pathogenic variants identified in the current study and in literature. Families A and B were described by Fattahi et al.[5]; individual UPN-0706 was described by Monies et al.[6] aa = amino acid. (B) Homology model of Znf domain 11 and 12 (in grey) in the human ZTBT11 in complex with DNA. Thr890 and Arg912 (in red) are located in the β-hairpin of Znf domain 11 and in the linker region between the domains, respectively.

The molar tooth-like sign seen on brain MRI in Patient FA.II-3 was suggestive of Joubert syndrome, a primary ciliopathy. However, reanalysis the WES data did not identify any putative pathogenic variants in known or candidate ciliopathy genes. Also, we did not identify dysfunction in genes related to cilia and ciliopathies in RNA-seq and ChIP-seq data of Patient FA.II-3.

We then queried the GeneMatcher database,[30] identifying additional patients carrying ZBTB11 variants (Families B and C). Exome sequencing in Family B identified compound heterozygous variants in ZBTB11 (GenBank: NM_014415.3) in the proband: Chr3(GRCh37):g.101384524T>A; c.907A>T; p.Ile303Phe maternally inherited and Chr3(GRCh37):g.101370393G>A; c.2779C>T; p.Arg927Ter paternally inherited with a CADD phred score of 26.9 and 37, respectively. The p.Ile303Phe was located in the BTB domain, while the p.Arg927Ter mapped in the Znf domain 12 (Fig. 2A). The ZBTB11 c.907A > T was not present in gnomAD, while c.2779C>T (rs767988210) was observed with allele frequency of 7.954e-6 and only in heterozygosity. Patients from Family C carried a homozygous missense variant in ZBTB11 (GenBank: NM_014415.3): Chr3 (GRCh37):g.101395605G>A; c.154C>T; p.Arg52Trp with a CADD phred score of 26.7 (Fig. 1S). This variant was not present in gnomAD and the amino acid change was located outside known functional domains (Fig. 2A).

We considered the ZBTB11 variants identified as putative disease causing, because they are rare in the population, are in silico predicted pathogenic, affect evolutionary conserved amino acid positions and segregated with the disease in the families. Also, bi-allelic ZBTB11 missense variants have been linked to intellectual developmental disorder, autosomal recessive 69 (MRT69) (OMIM 618383) in two families,[5] whose neurological presentation overlapped with the one of the patients in the three families. Another patient (ID: UPN-0706) has been reported with a homozygous ZBTB11 p.Arg912Gln and presenting global developmental delay, microcephaly, cataract, mild cerebellar atrophy and elevated methylmalonic acid.[6] This patient was described briefly as part of a large study, and additional clinical information, for example age or disease progression, was not available. Furthermore, knockdown of the fruit fly ZBTB11 orthologue CkIIα-i1 resulted in a neurological phenotype with reduced mushroom body, a neuropil structure of the brain.[5] Also, zebrafish expressing the mutant Zbtb11 p.C116S showed impaired craniofacial development and hydrocephalus with enlarged fourth ventricle,[17] resembling some of the findings in the patients.

Structure modelling of variants in ZBTB11

We performed structural analysis only for the variants located in the C-terminal Znf domain cluster, because no homology models were available for the analysis of the Arg52Trp and Ile303Phe variants. A structural model for the last two Znf domains (Znf domains 11 and 12) in the human ZTBT11 was generated from the experimental X-ray structure of the multidomain CCCTC-binding factor (CTCF) in complex with a known CTCF-binding site (PDB ID: 5UND)[35] using standard comparative modelling with SWISS-MODEL.[36] A model of the two Znf domains in complex with DNA was obtained by superposition with the CTCF(ZnF 4-10)-28-mer DNA complex. From the homology model (Fig. 2B), we can see that Thr890 is located within the β-hairpin of Znf domain 11 and Arg912 is located within the linker region between Znf domains 11 and 12. Linkers do play an active role in DNA binding.[37] The model suggests that Arg912 is interacting with the DNA backbone. Furthermore, the variant Chr3(GRCh37):g.101370393G>A; c.2779C>T results in a premature stop codon (p.Arg927Ter) within Znf domain 12, which will disrupt the folding of this domain if the transcript evades non-sense mediated mRNA decay. It is thus likely that two of the missense variants identified, p.Thr890Ala and p.Arg912Trp, and the premature stop codon (p.Arg927Ter) will lead to loss of function.

ZBTB11 binds genes encoding proteins involved in mitochondrial functions and in non-coding RNA processing

To identify genomic regions targeted by ZBTB11, we performed ChIP using an antibody against ZBTB11 on skin derived fibroblasts from six controls followed by sequencing (ChIP-seq: metrics available in Supplementary Table 1). ChIP-seq data from the six controls were analysed to identify consensus peaks, defined as overlapping peaks detected in at least four out of six samples. We identified 465 consensus peaks. When excluding peaks in unmapped genomic regions and in genes without identifier from Human Genome Organisation for Gene Nomenclature Committee (HGNC), the remaining 295 peaks were located within 238 genes (Supplementary Table 2). The 295 peaks mapped in promoters (82.4%), 1st exon (0.3%), 1st intron (3.7%), other introns (3.7%) and distal intergenic regions (9.8%) (Supplementary Fig. 1).

Gene ontology enrichment analysis of biological process terms connected to the 238 genes bound by ZBTB11 showed over-representation of biological processes related to protein translation, in particular mitochondrial translation and to non-coding RNA processing (Fig. 3 and Supplementary Table 3). Over-represented biological processes were related to the translation (GO:0006415; GO:0006414), including ribosome biogenesis (GO:0042254) and mitochondrial translation (GO:0006415; GO:0006414; GO:0070125; GO:0070126; GO:0032543; GO:0140053). Another cluster of over-represented biological processes was associated with processing of mRNA and non-coding RNAs, for example tRNA modifications and processing (GO:0002097, GO:0006400, GO:0008033, GO:0006399; GO:0006388); rRNA processing (GO:0006364; GO:0016072), RNA splicing (GO:0000394) and other non-coding RNA processing (GO:0034470). We also identified genes related to DNA conformation change and DNA duplex unwinding (GO:0032508; GO:0032392; GO:0071103), therefore related to transcriptional regulation.

Figure 3

GO enrichment analysis of biological process terms among the genes targeted by ZBTB11 in ChIP-seq experiment in control fibroblasts. ChIP-seq data from controls showed that ZBTB11 targets genes encoding proteins involved in mitochondrial functions and in RNA processes. Gene concept network plot (cnet plot) of the 238 genes associated to the ChIP-seq consensus peaks in controls depicts the linkages between genes and biological processes. Two main clusters of biological processes were identified: mitochondrial gene expression and RNA metabolic processes. The node sizes are dependent on the number of genes associated with the biological process as indicated. We used ClusterProfiler (R Bioconductor package) enrichment tool to perform analysis of over-representation of GO terms linked with our list of 238 significant genes identified in association with the differentially bound peaks. The GO enrichment analysis of biological process terms was subsequently visualized using the gene concept network plot (cnet plot).

Mutated ZBTB11 showed reduced binding to target genes

To assess the functionality of the mutated ZBTB11, we performed ChIP-seq in fibroblasts from Patients FA.II-1, FA.II-3 and FB.II-1 and compared with ChIP-seq data from the six control fibroblasts. We identified 70 differentially bound genomic peaks (P < 0.05), located in 61 genes, 29 of those present among the 238 genes identified by the consensus peaks (Supplementary Table 4). Forty-nine peaks were within promoters, eight within introns and four were in distal intergenic regions. The ChIP-seq data showed lower read counts at all 70 differentially bound genomic sites in the patients compared to the six controls (Table 2), but was not accounted for by differences in coverage (Supplementary Table 1, ChIP-seq quality metrics). The lower read counts indicated reduced binding ability of the mutated ZBTB11. Hence, the ZBTB11 variants in the patients behaved as loss of function.

Table 2

Total read counts of the significantly differentially bound genomic sites in the ChIP–seq data

Gene ID	Read depth controls	Read depth patients	Gene ID	Read depth controls	Read depth patients
ABHD13	2606.00	189.55	MPHOSPH10	1228.44	22.17
AC027119.1	169.09	0.00	MRPS34	6265.04	592.97
ACSF3	10171.97	54.26	MTERF4	4934.29	146.33
AGK	4403.60	194.95	MTG1	1217.51	18.27
AL162590.1	169.02	0.00	MTHFD2L	979.74	1.74
ANKRD26P1	263.24	0.00	NDUFS7	9178.26	847.26
ANTKMT	2214.36	118.60	OR7E116P	311.82	0.00
ATXN7	968.02	12.61	PCID2	8934.37	537.40
AURKAIP1	6134.92	222.58	PHLPP1	258.13	0.00
C8orf33	1179.67	19.23	PHRF1	222.53	0.00
C8orf33	5031.07	423.07	PIBF1	224.95	0.00
CLOCK	1043.13	29.70	PRKAR2A	213.82	0.00
CSMD3	715.05	0.00	RAD54L2	5400.04	211.47
CYREN	4051.38	296.40	RNF38	187.51	0.00
DDX11L8	744.13	0.00	RP11-618K16.4	1816.91	66.99
EEFSEC	3732.28	26.52	RP11-864G5.3 [Clone-based (Vega) gene]	147.64	0.00
ERI1	1346.47	44.90	RPL34P34	199.68	0.00
FANCA	903.03	0.00	SAFB	331.39	0.00
FBXW2	754.00	0.00	SDHAF3	1768.05	20.43
GLI3	231.76	0.00	SEC13	4427.94	138.51
GTF3C5	3184.77	293.99	SHARPIN	8191.02	717.46
GTPBP3	2625.46	263.40	SLC20A2	1385.28	56.47
HMGXB3	987.21	13.47	SLC25A53	776.27	15.65
JPX	3403.28	257.09	STX18	6114.85	646.28
KCTD5	3711.98	470.11	TEDC1	2204.63	89.21
LHFPL4	199.17	0.00	TERF2	142.59	0.00
lncRNA	306.63	0.00	TMEM129	1830.11	30.68
LRRC2	354.83	0.00	TMEM216	264.72	0.00
LRRC8A	1210.75	33.67	TMEM41A	3891.63	151.66
MAP2K5	891.70	4.09	TRDMT1	3594.15	310.06
MARK2P10	188.38	0.00	TRMU	4087.91	366.37
MCAT	3128.29	157.40	VWA8	2938.58	29.35
MCM3AP	2987.96	50.86	YY1	1281.07	34.72
MHENCR	1018.21	0.00	ZMYM4-AS1	206.55	0.00
MIR663AHG	341.73	0.00	ZSCAN22	3084.00	183.60

Distribution of the read counts in ChIP-seq data for the significantly differentially bound genomic sites (P < 0.05) between the six controls and the three patients (Patients FA.II-1, FA.II-3 and FB.II-1). The total read count was higher in the control group compared to the patient group for all differentially bound sites. The read numbers were obtained using the DiffBind (R Bioconductor package) with dba.peakset function. The values returned are controlled by the score parameter passed to dba.count. The default values are trimmed mean of M values (TMM) normalized, taking library size into account.

Gene ontology enrichment analysis was performed to identify of biological process terms connected to the genes harbouring the differentially bound genomic sites in the patient cells. This analysis revealed decreased binding of the mutated ZBTB11 to genes encoding for proteins involved in mitochondrial gene expression, and proteins related to RNA processing, modification and transport, including modification of tRNAs (Table 3). In addition, other genes showing reduced binding by the mutated ZBTB11 encoded mitochondrial enzymes (Table 2): acyl-CoA synthetase family member 3 (ACSF3), methylenetetrahydrofolate dehydrogenase activity (MTHFD2L), acylglycerol kinase (AGK), malonyl-CoA:ACP acyltransferase (MCAT), succinate dehydrogenase assembly factor 3 (complex II) (SDHAF3) and the NADH dehydrogenase (ubiquinone) Fe-S protein 7 (NADH-coenzyme Q reductase) (NDUFS7).

Table 3

Enriched biological processes of the genes showing decreased read counts in ChIP-seq data in cells from patients compared to controls

GO ID	Description	Gene ratio	P-value	P adj.	Gene ID
GO:0043648	Dicarboxylic acid metabolic process	3/52	0.002	0.186	ACSF3/SDHAF3/MTHFD2L
GO:0140053	Mitochondrial gene expression	4/52	0.001	0.186	AURKAIP1/MTG1/MRPS34/MTERF4
GO:0023035	CD40 signalling pathway	2/52	0.0007	0.186	FANCA/SHARPIN
GO:0016973	Poly(A)+ mRNA export from nucleus	2/52	0.001	0.186	PCID2/MCM3AP
GO:0050657	Nucleic acid transport	4/52	0.002	0.186	TERF2/PCID2/MCM3AP/SEC13
GO:0050658	RNA transport	4/52	0.002	0.186	TERF2/PCID2/MCM3AP/SEC13
GO:0051236	Establishment of RNA localization	4/52	0.002	0.186	TERF2/PCID2/MCM3AP/SEC13
GO:0006403	RNA localization	5/52	0.0004	0.186	TERF2/PCID2/MCM3AP/YY1/SEC13
GO:0034470	ncRNA processing	6/52	0.0007	0.186	TRDMT1/ERI1/MTERF4/TRMU/MPHOSPH10/GTPBP3
GO:0034660	ncRNA metabolic process	6/52	0.002	0.186	TRDMT1/ERI1/MTERF4/TRMU/MPHOSPH10/GTPBP3
GO:0001510	RNA methylation	3/52	0.001	0.186	TRDMT1/MTERF4/GTPBP3
GO:0009451	RNA modification	4/52	0.001	0.186	TRDMT1/MTERF4/TRMU/GTPBP3
GO:0006400	tRNA modification	3/52	0.001	0.186	TRDMT1/TRMU/GTPBP3
GO:0002098	tRNA wobble uridine modification	2/52	0.001	0.186	TRMU/GTPBP3
GO:0002097	tRNA wobble base modification	2/52	0.001	0.186	TRMU/GTPBP3

Enriched biological processes of the genes showing reduced binding in ChIP-seq data analysis in cells from Patients FA.II-1, FA.II-3 and FB.II-1 compared to six controls.

Transcriptome analysis revealed differential expression of genes related to mitochondria and CNS functions

To assess the extent to which the observed reduction in DNA binding ability of the mutated ZBTB11 resulted in transcriptional dysregulation, we performed RNA-seq of the fibroblasts from Patients FA.II-1, FA.II-3 and FB.II-1 and the six controls (RNA-seq metrics available in Supplementary Table 1). Analysis of the RNA-seq data identified 81 DEGs in the patients compared to the controls (P < 0.05), 51 were downregulated and 30 were upregulated (Supplementary Table 5). Among the 81 DEGs, EEFSEC was identified among the differentially bound peaks in ChIP-seq data in the patients, while EEFSEC, PIH1D2, TACO1, THEM4 and ZNF165 were also identified among the consensus peaks in ChIP-seq data in the controls. The DEGs identified were further analysed in GSEA[16] to identify dysregulated biological functions, which were found to be mostly related either to mitochondrial function or to neurons, dendrites and synapses (Fig. 4 and Supplementary Table 6). For example, the following DEGs were related to CNS (Fig. 4): CACNG8, MAP2, RNF157, VPS51. The Calcium Voltage-Gated Channel Auxiliary Subunit Gamma 8 (CACNG8) regulates trafficking of AMPA-selective glutamate receptors to cell membrane and synapses and modulates their gating properties.[38] CACNG8 is expressed primarily in hippocampus, cortex and subcortical regions, where it contributes to regulation of the social behaviour.[39] The microtubule associated protein 2 (MAP2) contributes to dendritic formation and maintenance, and in neural polarization and migration.[40] RNF157 is an E3 ubiquitin ligase expressed almost exclusively in the brain, shown to prevent apoptosis and promote survival of neurons, also required for dendrite growth and maintenance.[41] The Unc-13 homologue A (UNC13A) is involved in neurotransmitter release by acting in synaptic vesicle priming before vesicle fusion.[42]VPS51, downregulated in the patients, encodes for a vacuolar protein involved in retrograde transport from endosomes to the trans-Golgi. VPS51 loss of function causes pontocerebellar hypoplasia type 13 (PCH13; OMIM 618606),[43] which includes thinning of the corpus callosum, loss of white matter and cerebellar atrophy in common with MRT9.

Figure 4

GO enrichment analysis of biological processes and relative genes identified in the DEGs in RNA-seq data from patients and controls. Each horizontal line represents a GO term, while the vertical lines depict the genes detected in the RNA-seq experiments that clustered in each GO term. A gene may belong to multiple GO annotation categories. We used ClusterProfiler (R Bioconductor package) enrichment tool to perform over-representation analysis on GO terms linked to the DEGs in fibroblasts from patients versus controls. The GO enrichment analysis of biological process terms was subsequently visualized using the heatplot function. Only the DEGs included in the over-represented GO terms are shown.

Among the mitochondrial functions, the most prominent finding was the dysregulation of genes encoding core subunits of the mitochondrial respiratory complex I (ND2, ND3, ND4, ND5, ND6) and complex III (CYTB). In addition, several genes related to mitochondrial tRNAs were differentially expressed (Supplementary Tables 5 and 6), which conceivably impair translation: the selenocysteine-tRNA specific elongation factor (EEFSEC), the mitochondrially encoded tRNA-Glu (GAA/G) (MT-TE), the translational activator of cytochrome c oxidase 1 (TACO1), and the mitochondrial aminoacyl-tRNA synthetase (TARS2). We verified the reduced expression of EEFSEC and TACO1 using qPCR (Supplementary Fig. 2). Moreover, in line with the reduced binding of the mutant ZBTB11 to the ACSF3 promoter, we verified by qPCR decreased ASCF3 expression explaining the CMAMMA phenotype in patients (Supplementary Fig. 2). Both the ChIP-seq and RNA-seq data suggested a possible mitochondrial dysfunction in the patients. However, measurement of the OXPHOS enzyme activities in the patient fibroblasts did not reveal deficient enzyme activities compared to controls (Supplementary Table 7). During isolation of mitochondria from frozen samples of fibroblasts, part of the citrate synthase may leak out of the mitochondria and thus be lost. For this reason, we used complex IV as a mitochondrial marker enzyme. Normal OXPHOS enzyme activities in cultured fibroblasts do not exclude deficiencies in other energy-demanding tissues such as brain and muscle.

Discussion

We report five patients from three families with novel bi-allelic variants in ZBTB11. The variants met several pathogenicity criteria: in silico prediction of pathogenicity, segregation of the putative pathogenic variants with the affected individuals, phenotypic match among the affected individuals in the study and in literature,[5,6] and low allele frequencies compatible with a rare Mendelian disease. In addition, the ChIP-seq experiments showed reduced binding of the mutated transcriptional regulator ZBTB11 to target genes demonstrating a loss of function effect. We therefore concluded that the ZBTB11 variants detected were pathogenic. Dysfunction of ZBTB11 causes MRT69, a disease with autosomal recessive inheritance described in three families.[5,6]

Comparison of the clinical course of the five patients in the three families in this study with the patients in the two previously reported families with MRT69,[5,6] highlighted microcephaly, global developmental delay/intellectual disability, ataxia, hypotonia, spasticity, involvement of cranial nerves causing drooling and maxillary hyperplasia as shared features in this clinical entity. In addition, brain MRI showed a spectrum of abnormalities varying from pronounced cerebral cortex and cerebellar atrophy, ventricular dilatation and loss of almost all white matter seen in Family A, to isolated thinning of the corpus callosum seen in Family C. The published MRT69 patients exhibited cerebellar hypoplasia and mild ventriculomegaly, but without reduction of the central white matter.[5] In Family A, the disease was closely monitored over nearly four decades. Patients FA.II-1 and FA.II-3 appeared normal during the first months, but then started to deteriorate and by adolescence they were not able to walk and talk. Even though MRI examinations suggested loss of neural cells over the years, the neurological deficits in the two siblings remained rather unchanged after the regression during childhood. The rather stable clinical course observed between age 10 years and present day, suggested that most of the extensive atrophy likely occurred at an early age, which is in line with the MRI observations. It is reasonable to speculate that the atrophic process corresponded to the rapid deterioration of the clinical condition in late infancy, a feature overlapping with cases of mitochondrial-related leukodystrophy.[44] Another overlap with mitochondria-related disease were short stature and cataracts.[45,46] The patients in Families B and C also had progression of the disease, manifesting as language deterioration, although they retained the ability to speak and walk. Therefore, MRT69 is a progressive disease where clinical manifestations and severity vary.

The spectrum of disease’s severity suggest that different variants impair ZBTB11’s functionality to different degrees. In Family A, Thr890Ala is located in Znf 11 and Arg912Trp between Znf 11 and 12, in a linker that could play an active role in DNA binding.[37] Both variants could have a major impact on the functionality of ZBTB11, in line with the severe phenotype in the patients in Family A. Based on the limited clinical description available, Patient UPN-0706, homozygous for p.Arg912Gln located in the linker, also exhibited core manifestations of MRT69: global developmental delay, microcephaly, cataract, mild cerebellar atrophy and methylmalonic aciduria.[6]

While all other ZBTB11 variants identified in this study and in the literature[5,6] were missense and mapped in Znf domains or outside domains, Patient FB.II-1 had one stop-gain variant, Arg927Ter in Znf 12 and one missense, Ile303Phe in the BTB domain. In the RNA-seq data of Patient FB.II-1 we identified seven reads of the allele expressing Arg927Ter and 11 reads of the allele expressing Ile303Phe (data not shown), indicating that the allele expressing Arg927Ter to some degree escaped non-sense mediated decay and could encode a truncated protein with impaired folding of Znf 12. This may contribute to the less severe, but still progressive, disease in Patient FB.II-1. The BTB domain of a zinc finger protein can dimerize by binding the BTB domain of another zinc finger, modulating the ability of the zinc fingers to bind target DNA sequences.[47] Alternatively, a BTB domain can bind proteins without BTB domains, for example corepressors, regulating the transcriptional activity of the BTB-zinc finger protein.[47] The patients previously described with MRT69 and homozygous missense variants in Znf domains presented most of the pathological features so far documented in this syndrome: intellectual disability, microcephaly, delayed motor milestones, slight facial hypotonia and maxillary hyperplasia, spasticity, ataxia, cerebellar atrophy, thin corpus callosum and mild ventriculomegaly.[5] In Family C, the homozygous missense outside ZBTB11 functional domains might result in a mutated ZBTB11 with better preserved functionality. In line with the possibility, Family C had the mildest clinical presentation: mild intellectual disability, preserved speech (although with dysarthria), without spasticity, dystonia, only thinning of the corpus callosum and no CMAMMA. It should, however, be noted that the most affected patients, in Family A, were older than the patients in Families B and C, and further progression of their disease cannot be excluded.

The rare biochemical phenotype of CMAMMA was detected in the patients in Families A and B, but not in Family C, which had overall the mildest disease presentation. Elevated methylmalonic acid was reported in Patient UPN-0706,[6] while CMAMMA was not investigated in the patients reported by Fattahi et al.[5] (Dr Zohreh Fattahi, personal communication). In the patients, we could not identify putative pathogenic variants in ACSF3, however, we measured reduced binding of the mutated ZBTB11 to the ACSF3 promoter and reduced ACSF3 transcript levels. This probably caused reduced enzyme activity with subsequent accumulation of metabolites, manifesting as CMAMMA in the patients.

The molar tooth-like sign on brain MRI in Patient FA.II-3 was not seen in any of the other patients with ZBTB11 dysfunction, including those reported by Fattahi et al.,[5] whose brain MRI images were retrospectively analysed (Dr Zohreh Fattahi, personal communication). WES reanalysis did not reveal any likely pathogenic variants possibly explaining a ciliopathy in Patient FA.II-3, who did not present with polydactyly or other skeletal and visceral organ ciliopathy manifestations. Therefore, we cannot conclude whether this brain anomaly linked to Joubert syndrome was part of the MRT69 syndrome or not.

Both ChIP-seq and RNA-seq data analyses showed that ZBTB11 contributed to the regulation of mitochondrial functions, but the overlap in terms of genes identified in these two sets of experiments was limited, possible due to the low number of patients analysed expressing different ZBTB11 variants, and because of the cell type used. Also, fibroblasts from the patients showed OXPHOS enzyme activities comparable to those of the controls. Cell types more relevant to the clinical phenotypes under study might be more informative and should be used in additional functional studies addressing the mitochondria.

We documented that the mutated ZBTB11 showed reduced binding to target sites in the genome and that most of the DEGs were downregulated, confirming that ZBTB11 functions as a transcriptional activator and its loss of function primarily causes downregulation of gene expression.[48]

In controls, ChIP-seq data analysis showed that many genes targeted by ZBTB11 were involved in mitochondrial translation, and in processing of rRNA and tRNA, and other non-coding RNA species. In line with this, ChIP-seq data analysis of the patient cells showed a decreased binding of genes involved in expression of mitochondrial proteins and in processing of tRNA and other non-coding RNAs.

In the transcriptome analysis in patient cells, many DEGs encoded for mitochondrial proteins: for example, subunits of the Mitochondrial respiratory complex I (ND2-6) and III (CYTB), mitochondrial tRNAs (TARS2, MT-TE, EEFSEC). Another interesting observation in the RNA-seq data is the increased level of transcripts encoding for proteins involved in maturation and homeostasis of neurons, dendrites and synapses (UNC13A, MAP2, RNF157, CACNG8, VPS51), as dysregulation of these genes could have a role in the onset of the neurological phenotypes in the patients.

We documented that ZBTB11 has a role in ribosome biogenesis and rRNA processing, fitting with predominant nucleolar localization of ZBTB11 shown by Fattahi et al.,[5] suggesting its involvement in rRNA and ribosomal biogenesis processes. Overall, we showed reduced DNA binding by the mutated ZBTB11 and downregulation of genes involved in mitochondrial respiratory complex biogenesis, in mitochondrial translation and in processing of different RNA species in mitochondrial translation and in processing of different RNA species. All these findings are in line with the observations by Wilson et al.[48] in murine embryonic stem cells.

Conclusion

We describe five patients in three families with novel bi-allelic ZBTB11 variants presenting with clinical features compatible with MRT69. We expand our understanding of the varied phenotypic spectrum of MRT69. Our studies in patient fibroblasts documented reduced binding of the mutated transcriptional regulator ZBTB11, resulting in abnormal transcriptional regulation. Our data suggest affection of mitochondrial functions and RNA processing, contributing to the pathological phenotypes of the patients.

Click here for additional data file.