RRM1 variants cause a mitochondrial DNA maintenance disorder via impaired de novo nucleotide synthesis

Mitochondrial DNA (mtDNA) depletion/deletions syndromes (MDDS) encompass a clinically and etiologically heterogenous group of mitochondrial disorders caused by impaired mtDNA maintenance. Among the most frequent causes of MDDS are defects in nucleoside/nucleotide metabolism, which is critical for synthesis and homeostasis of the deoxynucleoside triphosphate (dNTP) substrates of mtDNA replication. A central enzyme for generating dNTPs is ribonucleotide reductase, a critical mediator of de novo nucleotide synthesis composed of catalytic RRM1 subunits in complex with RRM2 or p53R2. Here, we report 5 probands from 4 families who presented with ptosis and ophthalmoplegia as well as other clinical manifestations and multiple mtDNA deletions in muscle. We identified 3 RRM1 loss-of-function variants, including a dominant catalytic site variant (NP_001024.1: p.N427K) and 2 homozygous recessive variants at p.R381, which has evolutionarily conserved interactions with the specificity site. Atomistic molecular dynamics simulations indicate mechanisms by which RRM1 variants affect protein structure. Cultured primary skin fibroblasts of probands manifested mtDNA depletion under cycling conditions, indicating impaired de novo nucleotide synthesis. Fibroblasts also exhibited aberrant nucleoside diphosphate and dNTP pools and mtDNA ribonucleotide incorporation. Our data reveal that primary RRM1 deficiency and, by extension, impaired de novo nucleotide synthesis are causes of MDDS.


Multiple mtDNA deletions
The presence of multiple mtDNA deletions in patient muscle was assessed by Southern blot analysis. DNA was digested with a restriction enzyme that linearizes human mtDNA (PvuII), electrophoresed through agarose gels, and hybridized to a probe detecting mtDNA (Gene Link GeneProber GL557, #40-2055-41). Long range PCR protocols to detect multiple mtDNA deletions were previously described (1).

Whole exome sequencing
Family 1: WES was obtained from blood using Agilent SureSelectXT Human All Exon V5+UTRs capture and Illumina HiSeq2500 sequencing technology and analyzed for the presence of pathogenic mutations using the NextGENe software from Softgenetics and a proprietary analytical pipeline by Molecular Pathology of Columbia University.
Proband 3: WES was performed as previously described (2). Briefly, exome capture was attained using the Agilent Sure Select Human All Exon V5 (50 Mb) capture kit and sequenced using the Illumina GAIIx platform in 75 base pair reads. Reads were aligned to the human reference genome (UCSC hg19). Called variants were restricted to exonic (coding) or splice-site variants with a minor allele frequency of 0.01 (1%) or less from 378 in-house controls and external variant databases (GnomAD, NHLBI ESP, 1000 Genomes). Given the positive family history, autosomal dominant (heterozygous) inheritance was prioritized. Rare variants of nuclear genes associated with known mtDNA maintenance disorders were first analyzed. Next, variants were filtered using Gene Ontology (GO)-Terms to prioritize genes encoding products involved in mitochondrial function, DNA replication, or maintenance. GO-Terms employed were the wildcard term 'mitochond*', 'DNA repair', 'replication', 'transcription', 'nucleotide', 'purine', 'pyrimidine', 'exonuclease', 'polymerase', 'topoisomerase', 'ligase', 'helicase' and 'nucleoside'.
Copy number variants (CNVs) were also examined using the same GO-Terms. The identified RRM1 variant was confirmed by Sanger sequencing.

Magnetic resonance imaging
MRI of the brain of proband 1a without intravenous contrast was performed using multisequence multiplanar technique.
For MD simulations of the RRM1 variants, we used the X-ray structure of the human RRM1-TTP-GDP complex in its dimeric form at 3.21Å resolution (PDB ID 3HND). The RRM1 dimer consists of chains A and B, with residues 14-742 and 2-742, respectively, including Mg 2+ , GTP, TTP, and crystal water.
There are two cysteine residues, C352 and C356, close to the p.R381C variant. We modeled the C381-C356 disulfide bond for three reasons: (1) the sequence alignment between hRRM1 and ScRR1 in the Supplemental Figure 2A showed C352 (hRRM1) and S352 (ScRR1), and C356 (hRRM1) and A356 (ScRR1), and the hydrophobic property of Cys matched more closely with Ala than polar residue Ser; (2) both cysteine residues (C352 and C356) are in the loop region and have equal chances to form a disulfide bond with C381; and (3) in the WT there are hydrogen bonds between the guanidinium sidechain of R381 and the carboxyl group of E355, which makes R381 closer to neighboring C356. 4 The RRM1 was loaded into Maestro and processed using Schrödinger's Protein Preparation Wizard Panel (5) with default options such as missing hydrogen atoms and sidechain atoms and optimization of hydrogen bonding networks. Also, the optimal protonation states for histidine residues capped the N-and C-terminal ends with acetyl (ACE) and methyl amide (NMA) groups, respectively. Finally, the preprocessed structure was restraint minimized using an OPLS forcefield that minimizes hydrogen atoms while restraining the backbone and sidechain heavy atoms. The 3D Builder panel in Maestro was used to mutate residue R381 to C381 and H381. The final structure was further solvated with explicit water and ions in the cubic box for further molecular dynamics (MD) simulations.
We used all-atom explicit solvent MD simulations to investigate the conformational changes of RRM1 variants. The CHARMM-GUI solution builder (6) was used to prepare the system for these simulations. RRM1 was solvated in a cubic box with a dimension of 144.0 x 144.0 x 144.0 Å 3 . The system contained 0.15 M NaCl, and simulations were conducted at 310 K. The interactions between protein and water residues are modeled using CHARMM36 (7) force field and TIP3P (8) water. The forcefield parameters for TTP, GDP, and Mg 2+ were obtained from CGenFF (9). All simulations were run using GPU-version of GROMACS v5.1.4 (10). The simulation protocol started with energy minimization of 1000 steps, followed by constant atoms, volume, and temperature (NVT) equilibration (1 fs per time step) and NPT (2 fs per time step) simulations for 5 ns before start of the production simulations. During the initial minimization and equilibration, the backbone, sidechain, and dihedral restraints were applied, and in the production, the run protein was allowed to move freely without any restraint.
Temperature and pressure during the equilibration phase were maintained at 310 K and 1 bar, respectively, using Berendsen thermostat (11) and barostat. For the production run, a Nose-Hoover thermostat was used to control the temperature using a collision frequency of 1.0 ps -1 .
Lennard-Jones and electrostatic interactions were calculated explicitly within a cutoff of 1.2 nm, 5 and long-range electrostatic interactions were calculated by particle mesh Ewald summation.
The equilibration simulations were done on the NVT ensemble, and production run at constant atoms, pressure, and temperature (NPT) ensemble for 100ns. The total MD simulation time for all variants was 300 ns. The RRM1 WT, R381C, and R381H mutants contained ~281K atoms.
MD simulation trajectory analysis and visualization were done using Visual Molecular Dynamics (VMD 1.9.3) (12). Root mean square deviation (RMSD) was calculated using the RMSD trajectory tool in VMD. The distance between the center of mass (COM) of atoms was calculated using the VMD in-house Tcl scripts. The solvent-accessible surface area (SASA) and buried surface area (BSA) was calculated using the VMD measure sasa command with a probe of radius of 1.4 Å. All other time series figures were generated using Python scripts.

Ribonucleotide reductase activity
Measurement of RNR activity was assessed in actively cycling fibroblasts with a protocol adapted from Jong et al (13). Briefly, and to note protocol adaptations, fibroblasts were suspended 50 mM Tris pH 7.6, 20% glycerol, 15 mM MgCl 2, and 4 mM DTT and lysed by sonication. 500 μg protein was precipitated with 60% ammonium sulfate. Supernatant was passed through a G25 Sephadex column (Sigma, G25150). Reaction was performed at 25°C for This was followed by addition of 2.8 U NDP kinase and incubation at 25°C for 60 min. [2-14 C] CTP was measured using the DNA polymerase extension assay described by Martí, et al (14).
To measure mitochondrial dNTP pools, actively cycling fibroblasts were collected and analyzed as previously described by Martí et al (14).

Alkaline hydrolysis and Southern blot
DNA isolated from fibroblasts was hydrolyzed with 0.3 M KOH for 4 h at 55°C. DNA was then run out on a denaturing gel (50 mM NaOH, 1 mM EDTA), followed by Southern blot detection of mtDNA (Gene Link GeneProber GL557, #40-2055-41).

Statistics
We used two-way ANOVA with one-subject per cell method to compare the mean outcome of each proband to the controls. The validity for such an approach relied on the plausible assumption that the group difference in mean of the outcomes of interest does not depend on experiment. Therefore, the model only included group and experiment as the main effect but no group by experiment interaction (as such an interaction cannot be determined from the study data even if it existed). Findings were considered as statistically significant if the corresponding p-values were less than 0.05. Statistical analysis was performed using SPSS Version 28.