Imbalanced mitochondrial function provokes heterotaxy via aberrant ciliogenesis

Cilia are tubulin-based extensions, which can be motile or immotile. Immotile cilia, also known as primary cilia, exist as monocilia emanating from most vertebrate cells upon entry into G0 (1). They possess sensory functions and serve as an important hub for signaling processes such as the Hedgehog (Hh) pathway (2). Activity of such signaling pathways is crucial during development (2) and is also needed for heart formation (3, 4). Motile cilia, on the other hand, are typically found in multitudes, for instance on cells lining the surface of trachea, oviduct, and brain ventricles, and mainly serve to move fluids (1). A particular subset of motile cilia are monocilia within the left-right organizer (LRO), a transient organ of laterality — the node in mice or the Kupffer’s vesicle (KV) in zebrafish. These LRO cilia beat in a coordinated fashion and generate a counter-clockwise fluid flow that is translated into a biochemical signal. It provides the basis for the establishment of left-right (L-R) asymmetry and thus the correct placement of visceral organs, known as situs solitus (5). Because of the different kinds of cilia as well as their divergent functions in tissues and organs, ciliopathies are a heterogeneous group of diseases. Interestingly, though, situs anomalies such as those that occur in heterotaxy are a feature shared by many ciliopathies (6). Heterotaxy is recognized as a cause for heart malformations and consequently congenital heart disease (CHD) (7). This has recently been corroborated in an ENU-screening approach that revealed that mutations in genes affecting the generation and function of cilia lead to CHD (8).

We and others previously found that the mechanistic mTOR pathway is involved in the establishment of L-R asymmetry (9–11). When the mTOR pathway is deregulated, nodal cilia length is altered and symmetry breaking is not established properly (9, 11). Biogenesis of mitochondria via PGC1A is one of the many processes regulated by the mTOR pathway (12). Furthermore, mitochondria are often spatially located proximal to the base of cilia, e.g. as observed in photoreceptors, tracheal epithelium of vertebrates, or ciliates such as Tetrahymena thermophila (13–15). Together, these factors prompted us to investigate whether mitochondria influence the formation of cilia and consequently cilia-dependent processes such as L-R asymmetry development.

We wondered whether mitochondria could be found spatially proximal to the base of nonmotile primary cilia as it has been observed for motile cilia or photoreceptors (13–15). We therefore analyzed human fibroblasts via cryo-electron tomography, which revealed similar colocalization of mitochondria and ciliary base (Supplemental Video 1; supplemental material available online with this article; https://doi.org/10.1172/JCI98890DS1). Moreover, we found that microtubules project from the mother centriole at the cilium base toward mitochondria, suggesting a communication between the 2 organelles (Supplemental Video 1).

We then assessed the content of mitochondrial DNA (mtDNA) in samples of a heterogeneous cohort of heterotaxy patients. We detected a marked decrease of mtDNA levels in these patient samples as compared with healthy subjects (Figure 1A and Supplemental Figure 1). Yet samples derived from patients suffering from other isolated CHDs without visceral situs anomalies or additional cilia dysfunction phenotypes showed mtDNA levels comparable to healthy subjects (Figure 1A and Supplemental Table 1). This suggests reduced abundance of mitochondria specifically in heterotaxy patients. We therefore sought to functionally assess the influence of mitochondria on asymmetry development. To this end, we pretreated developing zebrafish embryos with the isoflavone 3-(2,4-dichlorophenyl)-7-hydroxy-4H-chromen-4-one (DCHC) to boost mitochondrial biogenesis (16) and generate elevated mitochondrial levels during developmental stages when cilia are involved in symmetry breaking (outlined in Figure 1B). Alternatively, we used rotenone to inhibit mitochondrial respiration during these developmental stages (Figure 1B). For both approaches we verified the efficiency of the treatment (Supplemental Figure 2). Enhancement as well as reduction of mitochondrial function caused phenotypes such as pericardial edema and curved tail, which can also be observed in zebrafish models of cilia dysfunction (Figure 1, C–E). Furthermore, we detected situs anomalies resembling heterotaxy after both treatments when we analyzed heart looping (Figure 1, F and G) and pancreas placement (Supplemental Figure 3). Spatially restricted expression of leftward genes such as southpaw (spaw) precedes heart looping and organ placement (17). We therefore examined spaw expression in zebrafish embryos after DCHC or rotenone treatment and observed increased occurrence of spatially aberrant expression (Figure 1, H and I). To verify this finding, we also injected zebrafish embryos with either capped RNA encoding murine Pgc1a to induce mitochondria biogenesis or a previously described translation-blocking morpholino (MO) targeting Cox5aa (18) to reduce mitochondrial respiration. Both treatments phenocopied the drug treatments with respect to disturbed heart looping and pancreas placement (Figure 1G and Supplemental Figure 3) and spatially aberrant spaw expression (Figure 1I). We thus conclude that mitochondrial homeostasis is important for correct situs development in zebrafish.

Figure 1

Mitochondrial content and function are associated with heterotaxy and laterality defects. (A) Blood samples from heterotaxy patients (n = 69) contain less mtDNA than healthy controls (n = 78) or samples from patients without heterotaxy-spectrum CHDs (n = 20). Kruskal-Wallis test with Dunn’s multiple comparison; ****P < 0.0001 for both comparisons. (B) Pharmacological manipulation of zebrafish embryos to stimulate (DCHC) or inhibit (rotenone) mitochondrial function. (C) Representative images of DCHC- and control-treated embryos 52 hpf. Arrow indicates pericardial edema. Scale bar: 200 µm. (D) Altered mitochondrial function induces pericardial edema. Unpaired, 2-tailed Welch’s t test; DMSO versus rotenone (rot), N = 8, n = 184/184, *P = 0.0059; DMSO versus DCHC, N = 5, n = 157/162, **P = 0.0037. (E) Altered mitochondrial function induces curved tail. Two-tailed Mann-Whitney test; DMSO versus rotenone, N = 6, n = 141/157, *P = 0.0043; DMSO versus DCHC, N = 4, n = 137/143, **P = 0.0286. (F, G) Detection of cml2c mRNA in 52 hpf zebrafish embryos to visualize heart looping. V, ventricle; A, atrium. Correct D-loop, incorrect L-loop, and no looping are shown. One-tailed Fisher’s exact test: DMSO versus rotenone, N = 5, n = 89/103, **P = 0.0023; CTRL MO versus Cox5aa MO, N = 3, n = 131/99, ***P = 0.0001; DMSO versus DCHC, N = 4, n = 111/115, ****P < 0.0001; Gfp mRNA versus mPgc1a mRNA, N = 3, n = 72/79, *P = 0.0031. (H, I) Ambiguous spaw expression upon mitochondrial manipulation. Normal expression, left (cartoon); aberrant expression, right or both sides. One-tailed Fisher’s exact test: DMSO versus rotenone, N = 8, n = 117/137, ****P < 0.0001; CTRL MO versus Cox5aa MO, N = 4, n = 125/128, **P = 0.0047; DMSO versus DCHC, N = 3, n = 150/124, ****P < 0.0001; Gfp mRNA versus mPgc1a mRNA, N = 5, n = 111/101, ***P = 0.0001. In D and E, red line indicates median. In F and H, scale bar: 50 µm.

LRO cilia are instrumental for proper establishment of L-R asymmetry (19). We hence visualized cilia in the KV, the LRO of zebrafish, after manipulation of mitochondria. We boosted mitochondria by treatment with DCHC before the KV formed so that at the onset of KV formation and cilium formation mitochondrial function was elevated (Supplemental Figure 2, B and C). Alternatively, we performed KV-specific overexpression of murine Pgc1a. Both approaches reduced cilia length in the KV, with the pharmacological intervention also reducing the number of cilia per KV (Figure 2, A–D). DCHC treatment during somitogenesis led to the same defects as treatment during gastrulation stages (Supplemental Figure 4). In contrast, inhibition of respiration via rotenone treatment or knock down (KD) of Cox5aa did not change the number of cilia per KV, yet resulted in elongated KV cilia (Figure 2, E–H). Manipulation of mitochondria thus inversely affects the length of motile cilia in the KV. Length changes of nodal cilia can affect ciliary beating and also the establishment of appropriate fluid flow in the KV, which ultimately triggers aberrant symmetry breaking (9–11, 20–24). We therefore analyzed ciliary beating and flow patterns in the KV of zebrafish embryos (Figure 2, I–M). Determination of the ciliary beat frequency (CBF) revealed no significant changes in DCHC-treated embryos, only a tendency toward a slightly increased CBF (Figure 2L). Rotenone treatment, however, reduced CBF significantly (Figure 2M). Unexpectedly, in the light of shorter cilia, an increased flow velocity could be detected in the anterior area of the KV in DCHC-treated embryos and in the left KV area for DCHC and rotenone, suggesting a disturbed flow pattern in this area (Figure 2, I–K). The left side of the KV is particularly important for symmetry breaking, which colocalizes with the downregulation of dand5 mRNA (previously known as charon, a negative regulator of the nodal cascade [ref. 25]) specifically at the left edge of the KV (21, 23). Consistent with the dampening of CBF and changes in nodal flow, we observed an increased frequency of embryos with aberrant dand5 expression, when treated with DCHC or rotenone (Figure 2N). We conclude that altered mitochondrial function affects nodal cilia length and beating, which as a consequence disturbs fluid flow in the KV and early steps in symmetry breaking such as spatial expression of dand5. To test whether this is also true for other motile cilia, we turned to another model system. The ciliate Tetrahymena thermophila requires functioning motile cilia for a number of processes, including swimming or phagocytosis. Rotenone and DCHC treatment affected both processes in Tetrahymena (Supplemental Figure 5). Similar to the increased flow speed in the KV and potentially due to increased ATP availability, DCHC treatment increased both the number of moving cells as well as speed of swimming, while rotenone treatment reduced both these parameters. Perhaps this reflects cilia movement in the KV of zebrafish embryos treated with these substances (Figure 2, L and M). We conclude that formation and function of motile cilia depends on mitochondrial homeostasis and that this connection is conserved.

Figure 2

Dosage of mitochondrial function reciprocally influences length of LRO cilia and disturbs flow. (A–D) Increased mitochondrial biogenesis decreases LRO cilia length. acTub, acetylated tubulin, cilia. (C) Cilia per KV. DMSO versus DCHC, 19/16 embryos, ***P = 0.0002; Gfp versus mPgc1a mRNA, 7/13 embryos. (D) LRO cilia length. DMSO versus DCHC, N = 19/16, n = 768/339; Gfp versus mPgc1a mRNA, N = 7/13, n = 230/357. ****P < 0.0001. (E–H) Inhibition of mitochondrial respiration elongates LRO cilia. (G) Cilia number per KV. DMSO versus rotenone, 11/15 embryos; CTRL MO versus Cox5aa MO, 20/20 embryos. (H) LRO cilia length. DMSO versus rotenone, N = 11/15, n = 643/700; CTRL MO versus Cox5aa MO, N = 20/20, n = 529/473. ****P < 0.0001. (I) Heat maps of flow speed. Orientation and color code are indicated. (J) LRO-flow speed after DMSO or DCHC (N = 13 and 12). Overall n = 722/607, anterior n = 399/319, *P = 0.0108, posterior n = 323/288, left n = 381/314, **P = 0.0169, right n = 342/293. (K) LRO-flow speed after DMSO or rotenone (N = 9 and 8). Overall n = 353/400, anterior n = 177/186, posterior n = 176/214, left n = 175/187, *P = 0.0422, right n = 178/214. (L) Ciliary beat frequency (CBF) after mitochondria stimulation. DMSO, 13 embryos, 113 cilia; DCHC, 12 embryos, 109 cilia. (M) CBF after inhibition of mitochondria. DMSO, 8 embryos, 66 cilia; rotenone, 9 embryos, 73 cilia; ****P < 0.0001. (N) dand5 (charon) expression at 10 ss. Shown are correct expression (right) and aberrant expression (both, left). Two-sided Fisher’s exact test: N = 3; DMSO versus DCHC (n = 67/51), ****P < 0.0001; DMSO versus rotenone (n = 82/77), **P = 0.0083. Two-tailed Welch’s t test in C and G (left); 2-tailed Mann-Whitney test in D and G (right), H, and J–M. Scale bars in A, B, E, F, and N: 20 μm. Shown are mean ± SEM in panels C and H. Red line indicates median.

Nonmotile primary cilia, however, are also involved in the process of symmetry breaking. These cilia lining the KV are thought to sense the flow and translate it into leftward gene expression (26, 27). Therefore, we analyzed whether immotile cilia such as primary cilia generation is also sensitive to altered mitochondrial function. We manipulated mitochondria in human fibroblasts with the same compounds as in zebrafish (Figure 3, A–H, and Supplemental Figure 6). Change of mitochondrial function in any direction did not impact the number of ciliated cells (Figure 3, C and D). However, DCHC treatment reduced cilia length, whereas rotenone caused an elongation of primary cilia, as we have observed with motile LRO cilia (Figure 3, E and F). Length changes of primary cilia also affect their function. To assess this, we measured the ability of cells to induce Hh signaling upon stimulation with Smoothened agonist (SAG) after the cells have been pretreated with either DCHC or rotenone. Both pretreatments triggered reduced capacity for Hh signaling induction (Figure 3, G and H). This indicates that the functionality of primary cilia is reduced in cells with a state of altered mitochondrial function.

Figure 3

Dosage of mitochondrial function reciprocally influences length and function of primary cilia in human fibroblasts. (A, B) Increased mitochondrial function causes shortening of primary cilia of human fibroblasts. Representative images of cells treated with DMSO or DCHC. Cilia are visualized by staining for acetylated Tubulin (acTub, green), while COX IV (red) stains mitochondria. (C, D) DCHC or rotenone treatment does not change the number of ciliated cells. Paired 2-tailed t test; DMSO versus DCHC, N = 3, n = 248/246, P = 0.3517; DMSO versus rotenone, N = 3, n = 253/249, P = 0.1361. Shown are mean ± SEM. (E, F) DCHC or rotenone treatment alters the length of cilia. Two-tailed Mann-Whitney test; DMSO versus DCHC, N = 3, n = 96/105, ****P < 0.0001; DMSO versus rotenone, N = 3, n = 106/103, **P = 0.0011. Red line indicates median. (G, H) DCHC or rotenone treatment reduces the cellular ability to induce the Hedgehog-target gene GLI1 upon stimulation with an agonist of the Hh pathway. Paired, 2-tailed t test; DMSO versus DCHC, n = 4, *P = 0.0155; DMSO versus rotenone, N = 3, **P = 0.0277. Shown are mean ± SEM. (I, J) Patient-derived human fibroblasts carrying mutations in MPV17 or NFU1 extend longer primary cilia as compared with control cells (WT). Acetylated Tubulin (acTUB, red) stains cilia, while γTubulin (γTUB, green) stains centrioles. Kruskal-Wallis test with Dunn’s multiple comparison test, N = 3, n = 104/95/107, WT versus MPV17, **P = 0.0028; WT versus NFU1, ***P = 0.0002. Red line indicated median. In A, B, and I, scale bars: 10 μm.

To further test the influence of mitochondria on ciliogenesis, we analyzed primary cilia on cells derived from patients carrying biallelic pathogenic variants of MPV17 and NFU1, respectively, which cause mitochondrial DNA depletion syndrome 6 (phenotype MIM 256810) and multiple mitochondrial dysfunctions syndrome (phenotype MIM 605711), respectively (Figure 3, I and J) (28, 29). Both patient-derived cell lines showed longer cilia as compared with wild-type control cells, thus reflecting the same phenotype as observed after inhibition of mitochondrial respiration (Figure 3, I and J). Formation of primary and motile cilia therefore respond the same to altered mitochondrial function. We hence propose that mitochondrial activity is precisely controlled to ensure correct cilium formation in general, and that increases and decreases in mitochondrial activity influence cilia length in inverse directions, which consequently reduces cilia functionality.

Our finding that mitochondrial dysfunction confers defects in cilium formation suggests that mutations in mitochondria-associated genes could cause situs anomalies such as heterotaxy. Ciliopathy phenotypes are rare in mitochondriopathy patients, but have been reported in certain cases (30–34). Nevertheless, there is phenotypic overlap in symptomatology. For example, patients with syndromic ciliopathies often possess phenotypes such as renal cysts or retinal degeneration that also occur in mitochondrial disorders (35). To further test our hypothesis, we performed whole-exome sequencing in a cohort consisting of 285 heterotaxy patients and carried out a focused analysis of 1174 genes listed in the Human MitoCarta2.0 (36) and the BCM-MitomeNGS (Baylor Genetics Laboratories) (Supplemental Table 2). We found either homozygous or X-linked hemizygous mutations in mitochondria-associated genes in 26 patients (9.12%; Table 1 and Supplemental Table 3). Six of these genes (MRPL44, MRPS25, MRPL38, NDUFA1, COQ9, TAZ) have previously been linked to known oxidative phosphorylation (OXPHOS) diseases (37). Overall, the genes are involved in diverse processes connected to mitochondria, among them synthesis of cardiolipin (TAZ), mitochondrial translation (MRPL44, MRPS25, MRPL38), or mitochondrial ubiquinone biosynthesis (COQ9). Some of the detected rare variants have not been listed in population databanks such as the Exome Aggregation Consortium (ExAC) browser (MTRR, CA5B, IDH3G, SLC25A5, SLC25A53, TIMM17B). In addition, the minor allele frequency of variants present in ExAC indicates that they are ultra-rare and present at a significantly higher frequency in heterotaxy cases than within the general population (Table 1 and Supplemental Table 3). Next, we tested a more complex genetic inheritance model by investigating variant burden in order to analyze a potential association with heterotaxy. We considered the number of rare, predicted damaging variants in our cohort and 1000 Genomes controls to establish the variant burden. Examining 1174 mitochondria-associated genes identified a significant increase in damaging variant burden in the overall heterotaxy cohort as compared with 1000 Genomes controls (Figure 4A). When stratified by ancestry, patients with European and African ancestries showed a significantly increased burden of rare damaging variants (Figure 4A). To test the specificity of this finding, a comparison was performed against exon-length–matched control gene lists including autism, cancer, and random disease-associated genes (Figure 4B and Supplemental Tables 2 and 4). A significantly higher burden of rare damaging variants was identified in mitochondria-associated genes in the heterotaxy cohort.

Figure 4

Increased burden of rare damaging variants in mitochondria-associated genes in patients with heterotaxy. (A) Variant burden analysis in mitochondria-associated genes in the heterotaxy cohort versus 1000 Genomes controls. Burden is presented as ratio of damaging/tolerated rare variants. Wilcoxon rank sum test, 2-sample test with normal approximation; overall n = 2503/284, ****P < 0.0001; European (EUR), n = 502/190, *P = 0.0248; African (AFR), n = 661/41, **P = 0.0109; Admixed American (AMR), n = 347/32, P = 0.1502; Asian, n = 993/21, P = 0.106. Arrows indicate extreme outliers which are included in analysis but not shown. (B) Variant burden analysis in mitochondria-associated genes compared with control genes in the heterotaxy cases. Nonparametric comparisons for each pair using Wilcoxon method; n = 283; mitochondria-associated genes versus autism genes, ****P < 0.0001; cancer genes, ****P < 0.0001; random disease genes A, **P = 0.0012; random disease genes B, ***P = 0.0015. Data sets are represented in outlier box plots indicating median and 25th and 75th quantile, respectively. Whiskers indicate first quartile –1.5* (interquartile range) and third quartile +1.5* (interquartile range), respectively. Dots indicate outliers, and arrows indicate extreme outliers which are included in analysis but not shown.

Table 1

Rare variants in mitochondria-related genes detected in heterotaxy patients

To evaluate these mitochondria-associated genes, homologs of the homozygous and hemizygous patient variants were analyzed in zebrafish and human fibroblasts. Affected genes are expressed during early somitogenesis in zebrafish embryos when symmetry breaking is taking place (Supplemental Figures 7 and 8). To functionally test the influence of the detected variants on ciliogenesis, we chose 6 candidate genes on the basis of predicted damage to the protein and expression in human fibroblasts (judged by proteinatlas.org; ref. 38) namely MTRR, CA5B, ACOX1, TAZ, OSBPL1A, and PHKA2. Then we targeted the candidate genes by siRNA-mediated KD in human fibroblasts and confirmed KD efficiency (Supplemental Figure 9). KD of each candidate gene led to reduced mitochondrial function as measured by a lower ATP content (Figure 5A). Consequently, and in line with our previous results, primary cilia were elongated after KD of candidate genes as compared with control cells (Figure 5, B and C), whereas the number of ciliated cells remained unaffected (Supplemental Figure 10), indicating that affecting mitochondrial function via multiple pathways shares the outcome of elongated primary cilia in cells.

Figure 5

Cilia length and function defects upon loss-of-function of heterotaxy patient-derived genes. (A) KD of candidate genes reduces cellular ATP. One-way ANOVA with Holm-Sidak comparison, n = 10, all compared with control siRNA (siCTRL): ****P < 0.0001 (TAZ, OSBPL1A, MTRR, PHKA2), ***P < 0.0001 (CA5B), **P = 0.0049 (ACOX1). Red line indicates median. (B, C) KD of candidate genes causes longer cilia. Cartoon illustrates staining for cilia (acTub, red), the basal body (BB, green), and the nucleus (blue). Scale bar: 10 μm. Kruskal-Wallis test with Dunn’s multiple comparison, N = 3, n = 95/99/96/96/98/94/100, all compared with siCTRL: ****P < 0.0001 (CA5B, TAZ, OSBPL1A, MTRR), **P = 0.0049 (ACOX1). Red line indicates median. (D) Representative images of 48 hpf zebrafish embryos, which were either noninjected (NI), control-injected (CTRL MO), or injected with MO targeting splice sites (spMO) in Mtrr, Phka2, and Osbpl1a, respectively. (E) Quantification of pericardial edema and body curvature in embryos. One-way ANOVA with Holm-Sidak comparison; N = 7/7/3/4/4, n = 196/178/52/102/109, all compared with CTRL MO. Pericardial edema: ****P < 0.0001 (Mtrr splMO, Phka2 splMO); *P = 0.0144 (Osbpl1a splMO); curved tail: ****P <0.0001 (Mtrr splMO), ***P =0.0008 (Osbpl1a splMO), **P = 0.0003 (Phka2 splMO). (F) KD of candidate genes increases frequency of aberrant heart looping. Two-sided Fisher’s exact test; N = 9/9/3/6/4; n = 227/200/50/143/104; all compared with CTRL MO: ****P <0.0001 (Mtrr splMO, Phka2 splMO), ***P = 0.0006 (Osbpl1a splMO). (G) KD of candidate genes does not alter the number of cilia per KV. One-way ANOVA with Holm-Sidak comparison, n = 15/12/13/15. (H) KD of candidate genes confers elongation of LRO cilia. Kruskal-Wallis test with Dunn’s multiple comparison test; n = 261/209/268/343; all compared with CTRL MO, ****P < 0.0001. Red line indicates median. Scale bars: 10 μm (B) or 500 μm (D). Shown are mean ± SEM in E and G.

We next tested whether KD of these candidate genes in zebrafish embryos is sufficient to influence ciliogenesis and cause ciliopathy-like phenotypes. We therefore designed splice-blocking morpholinos (splMOs) for candidate genes for which it was possible, namely MTRR, OSBPL1a, PHKA2, or TAZ, and verified the targeting efficiency of the splMOs (Supplemental Figure 11). Targeting of all tested candidate genes resulted in similar phenotypes as in the initial mitochondria manipulation experiments (Figure 5, D and E, and Figure 6, A and B). Morphant embryos also showed an increased occurrence of aberrantly looped hearts and misplaced pancreas as compared with control-injected embryos (Figure 5F, Figure 6C, and Supplemental Figures 12 and 13) suggesting that all tested genes are important for normal L-R asymmetry development. This hypothesis was further strengthened by cilia measurements in the KV. Cilia of morphant embryos were elongated, whereas the cilia number per KV was unchanged (Figure 5, G andH, and Figure 6D). Hence, KD of candidate genes identified in heterotaxy patients compromises elongation of primary as well as nodal cilia as judged by their increased length, ultimately leading to disturbance of cilia-dependent processes such as symmetry breaking.

Figure 6

Patient rare variants cannot rescue ciliopathy-like phenotype. (A) Representative images of 48 hpf zebrafish embryos, which were either noninjected (NI), control-injected (CTRL MO), or injected with MOs targeting Taz (Taz MOs), respectively. Optionally, capped RNA encoding wild-type hTAZ or hTAZ-R253Q were coinjected with TAZ MOs. Scale bar: 500 μm. (B) Quantification of pericardial edema and curved tail after Taz KD. N = 5, n = 119/108/119/124/109; pericardial edema, RM 1-way ANOVA with Sidak’s multiple comparison, CTRL MO versus Taz MOs, **P = 0.0064; Taz MOs versus Taz MOs + hTAZ, *P = 0.0328; MOs + hTAZ versus Taz MOs + R253Q, ***P = 0.0037; curved tail, Friedman test with Dunn’s comparison, CTRL MO versus Taz MOs, ****P = 0.0055. (C) Rescue of heart looping. N = 5; n = 119/94/116/118/112; 2-sided Fisher’s exact test; CTRL MO versus Taz MOs, ****P < 0.0001; Taz MOs versus Taz MOs + hTAZ, ***P = 0.0002; Taz MOs versus Taz MOs + R253Q, P = 0.0828; Taz MOs + hTAZ versus Taz MOs + R253Q, ****P < 0.0001. (D) Altered expression of TAZ does not affect number of cilia per KV. Kruskal-Wallis test with Dunn’s multiple comparison, n = 8/8/8/7. Coinjection of RNA encoding rare variant of hTAZ as detected in heterotaxy patients does not rescue LRO cilia length in TAZ morphant zebrafish. Kruskal-Wallis test with Dunn’s multiple comparison; N = 8/8/8/7, n = 340/295/315/326; CTRL MO versus Taz MOs, ****P < 0.0001; Taz MOs versus Taz MOs + hTAZ, ****P < 0.0001; Taz MOs versus Taz MOs + R253Q, P = 0.737; Taz MOs + hTAZ versus Taz MOs + R253Q, ****P < 0.0001. Shown are mean ± SEM and the median as red line, respectively, in B and D.

To further corroborate these observations we performed rescue experiments. We treated TAZ KD human fibroblasts with mitoTempo, a mitochondrion-specific scavenger of ROS that has been shown to restore mitochondrial function in this setting (39). This also worked in our hands as judged by increased cellular ATP levels in TAZ KD cells (Supplemental Figure 14). Consequently, mitoTempo treatment of TAZ KD cells caused shortening of primary cilia, while the number of ciliated cells remained the same (Supplemental Figure 14). Thus, we propose mitochondrial dysfunction after TAZ loss-of-function in human fibroblasts as the likely cause for elongation of primary cilia.

Next, we aimed to rescue ciliopathy-like phenotypes in zebrafish. We knocked down Taz and Mtrr and coinjected capped RNA encoding either the wild-type protein or the variant as detected in the heterotaxy patients. KD of both candidate genes caused pericardial edema or curved tail (Figure 6B and Supplemental Figure 15). Moreover, whereas injection of RNA encoding the wild-type protein was able to at least partially rescue these defects, RNA encoding the patient rare damaging variants failed to do so (hTAZ-R253Q or hMTRR-R552Q; Figure 6B and Supplemental Figure 15). Coinjection of RNA encoding wild-type TAZ further restored normal heart looping and pancreas positioning in Taz morphant embryos, whereas RNA encoding hTAZ-R253Q did not (Figure 6C and Supplemental Figure 13). Similarly, wild-type MTRR restored L-R asymmetry in zebrafish depleted of Mtrr, but the patient variant hMTRR-R552Q did not (Supplemental Figure 16). Finally, although the number of cilia per KV was unchanged, KD of Taz in zebrafish embryos triggered elongation of LRO cilia, which was fully rescued by coinjection of RNA encoding hTAZ (Figure 6D). In contrast, the patient variant hTAZ-R253Q was not able to rescue LRO cilia elongation, suggesting that this variant at least contributed to the development of heterotaxy in the respective patient (Figure 6D).

So far, dysfunctional Taz protein has mainly been associated with Barth syndrome, a rare X-linked syndrome characterized by early-onset skeletal and cardiac myopathies due to mitochondrial remodeling (40). In zebrafish, loss of Taz had been associated with a general impairment of cardiac function such as pericardial edema, which we also observed, and bradycardia (41). We therefore assessed the heart rate in Taz morphant embryos and counted a slight but significant reduction in heartbeats per minute (Supplemental Figure 17). Only coinjection of RNA encoding the wild-type hTAZ rescued this phenotype, while coinjection of RNA encoding the patient variant hTAZ-R253Q failed to do so (Supplemental Figure 17). Therefore, our results suggest that variants of the TAZ gene have the potential to impact the organism at time points earlier than previously described and on multiple levels.

Mitochondria are historically best known for their role in energy production. Here we show that mitochondrial function is also connected to cilium formation and function. We found a significant reduction in mtDNA selectively in a cohort of heterotaxy patients and consequently manipulated mitochondrial function in zebrafish embryos, human fibroblasts, and Tetrahymena thermophila. In response to such manipulation, we detected heterotaxy-like phenotypes and identified cilia length aberrations, which inversely correlate with diminished or increased mitochondrial function. KD of mitochondria-associated candidate genes identified in a second cohort of heterotaxy patients attenuated ATP abundance in human fibroblasts and conferred cilium elongation in fibroblasts and zebrafish, resulting in ciliopathy-like phenotypes including L-R asymmetry defects. Patient variants failed to rescue KD-induced phenotypes, indicating that the identified rare variants may be contributors or even causes for the development of heterotaxy. Yet, how could reduced ATP levels cause elongation of cilia? Inhibition of energy-dependent processes that reduce cilia length possibly explains the observed phenotype. A candidate factor could be the kinesin KIF19A, which locates to the tip of cilia and depolymerizes ciliary microtubules from the plus end in an ATP-dependent fashion (42). Knock out of Kif19a consequently causes longer cilia in different organ systems of the mouse (42). Lower abundance of ATP due to mitochondrial dysfunction could reduce the processivity of Kif19a at the tips of cilia and hence attenuate Kif19a’s length-restricting activity, leading to elongated cilia. Vice versa, elevated ATP levels potentially increase Kif19a activity, leading to shorter cilia.

However, mitochondria have additional important functions besides energy production. Indeed, the rare genetic variants we detected in exome-sequencing data of heterotaxy patients include genes affecting very diverse mitochondrial functions. These include among others heme biosynthesis (ALAS2), DNA repair (APEX2), mitochondrial cardiolipin production (TAZ), mitochondrial membrane transport (TIMM17B), lipid metabolism and ROS regulation (ACOX1), mitochondrial protein translation (MRPL38, MRPL44, MRPS25), ubiquinone synthesis (COQ9, ADCK4), methionine synthesis (MTRR), glycogen storage (PHKA1, PHKA2), lipid transfer (OSBPL1A), and mitochondrial HCO₃^– production (CA5B). Yet despite their diverse functions, KD of all investigated candidate genes reduced ATP content and simultaneously caused longer cilia in our functional tests. This finding was further complemented by our observation of longer cilia in patient-derived human fibroblasts with defects in mitochondrial function and maintenance, respectively. Therefore, disturbances of mitochondria on very diverse functional levels seem to funnel into the same phenotype of elongated cilia. The question that remains is whether reduced ATP levels are the common trigger for cilia elongation. This could be a possibility, since KD of diverse candidate genes led to attenuated ATP levels. Moreover, different mitochondrial processes not directly involved in respiration have been shown to impact ATP synthesis: fatty acid metabolism, for instance, is coupled to ROS generation and the respiratory chain. In mouse cardiomyocytes, lipid overload results in increased respiratory activity and hence ATP production. Extended exposure to excessive lipids, on the other hand, results in the production of ROS (43). Calcium overload in mitochondria also impacts ATP production through mitochondrial permeability transition pore opening and subsequent membrane depolarization (44). Mitochondrial protein translation is coupled to ATP synthesis, too, as loss of the eukaryotic translation initiation factor 2B subunit 5 results in a reduction of respiratory chain components on a protein level and subsequently a decrease in ATP synthesis (45). Similarly, a very recent, large-scale Crispr/Cas9 approach identified deficiencies in mitochondrial protein translation as one of the prime causes of impaired respiratory chain activity although the 5 complexes of the respiratory chain are largely encoded by nuclear genes (37). Interestingly, 3 genes for which we detected homozygous rare damaging genetic variants, namely MRPL38, MRPL44, and MRPS25, regulate mitochondrial translation and are known as OXPHOS disease genes (37). Methionine metabolism has also been linked to mitochondrial respiratory function. Methionine is enriched in the intramembraneous space and serves as antioxidant for the respiratory chain (46). The mitochondrial methionine cycles further provide the source of methyl groups for epigenetic changes in DNA (47). Hypermethylation of the Cox5a promoter has been linked to reduced expression of Cox5a, which we also depleted to inhibit respiratory activity in the LRO of zebrafish embryos, and thus reduced activity of the respiratory chain, resulting in less ATP synthesis (48). Hence, the assessment of cellular ATP content may serve as a readout of mitochondrial function, and potentially also as a diagnostic in afflictions connected to impaired cilia formation. Vice versa, manipulation of ATP content could possibly be a means of correcting cilia defects. Moreover, we measured strikingly attenuated mtDNA content specifically in peripheral blood samples of heterotaxy patients as opposed to samples of patients with CHD not connected to ciliary dysfunction or from healthy subjects. This might present as an additional diagnostic tool to recognize ciliopathies, as peripheral blood is easily accessible and measurement of mtDNA content straightforward. However, further studies are needed to scrutinize and verify the usefulness and robustness of such diagnostic and therapeutic approaches.

Cilia defects can give rise to CHD (8, 49). In addition, small studies and case reports have highlighted that certain cases of mitochondriopathies also display congenital heart malformations and patients with CHD may have reduced mitochondrial DNA content (50–52). Interestingly though, mitochondria and cilia organelles have not previously been associated mechanistically as a possible common link in CHD causation. For example, TAZ loss-of-function variants, which cause mitochondrial insufficiencies, have been associated with fetal and early postnatal onset cardiomyopathy, but not with cilia defects, although Taz knockout mice display spermatogenesis defects and Drosophila flies devoid of functional taz were sterile, both being defects suggesting potential abnormal cilia function (53, 54). In a recent ENU screen for novel predisposing genes for CHD (8), 12 of the 61 human homologs of the identified genes localize to mitochondria (according to genecards.org, ref. 55), which corresponds to 19.7% of the identified susceptibility genes. This suggests that pathogenic rare variants in genes connected to mitochondria can contribute to the development of CHD at least in the murine system. Moreover, at least one of the laterality genes (ATIC), for which genetic variations were reported in a large human cohort of CHD (49), also localizes to mitochondria.

A phenotypic overlap between mitochondrial disorders and cilia-related abnormalities exists (e.g., nephronophthisis, craniofacial anomalies, microcephaly, and obesity) (18, 30–35). A previous study has identified pathogenic variants of mtDNA at elevated frequencies in cardiac tissue from patients suffering from CHD (50). The connection of mitochondrial function and ciliogenesis identified by the current study raises the possibility that pathogenic variants in nuclear-encoded mitochondrial genes lead to reduced mitochondrial function and also defective cilia. These data together indicate that mutations in nuclear- as well as mitochondrial-encoded genes may contribute to defects triggered by dysfunctional cilia such as CHD. Thus, in order to detect the disease-causing rare variants in patients presenting with ciliopathy-like features, mitochondria-related genes should also be taken into consideration.

The data demonstrate that manipulation of mitochondrial function in multiple organisms has the potential to alter generation of cilia. Although syndromic and Mendelian forms of CHD explain a proportion of disease in humans, many cases of CHD behave as a complex trait. The increased burden of rare damaging variants in mitochondria-associated genes identified in our heterotaxy cohort suggests that the accumulation of protein-altering variants that collectively disrupt mitochondrial function act as potential susceptibility factors for heterotaxy and CHD.

Patient cohorts. Two heterotaxy cohorts were used in this study: To assess the content of mtDNA, DNA samples from heterotaxy patients (n = 69), patients suffering from non–heterotaxy-related CHDs (n = 20, Supplemental Table 1) and healthy subjects (n = 78, unaffected individuals, mostly parents of patients) provided by the German National Register for Congenital Heart Defects (Competence Network for Congenital Heart Defects, Berlin, Germany) were used. DNA samples were from whole blood and isolated by standard procedures. The local ethics committee in Berlin approved collection and distribution (due to German practice, no ethnicity data were collected). Medical records regarding the heart phenotype and situs were provided with the samples. The local ethics committee in Ulm approved mtDNA content analyses for this cohort.

Information from a second cohort (n = 285) with phenotype information including imaging and diagnostic studies, clinical genetic testing results, and pathology results was collected at Cincinnati Children’s Hospital Medical Center (CCHMC) under a protocol approved by the CCHMC Institutional Review Board. The cohort consisted of patients with a diagnosis of heterotaxy or heterotaxy-spectrum heart defects. DNA samples from whole blood were isolated by standard procedures.

Zebrafish husbandry and manipulation. Zebrafish wild-type strains AB and EK (gift of Kenneth D. Poss, Duke University Medical Center, Durham, NC) were maintained under a 14-hour light, 10-hour dark cycle and fed a standard diet. Embryos were generated by natural matings and allowed to develop to the desired stages at 28.5°C. Microinjections into zebrafish eggs were carried out with a Femtojet microcompressor (Eppendorf), pulled capillaries, and a Narishige micromanipulator. Antisense morpholino oligonucleotides (MOs) or capped RNAs were injected into the yolk at the 512–1000 cell stage to specifically target cells giving rise to the KV (56). The Cox5aa MO and the Taz ATG MO have been previously described (18, 41). Heterotaxy candidate genes were targeted with the following MOs: Mtrr splMO (5′-GCTTGAGAAAGATGTGGAAAAATAT), Osbpl1a splMO (5′-AAGGATGTGTTTTACCCTTGCAGTT), Phka2 splMO (5′-CCACCTGCAGTTACATTGACGTG), and Taz splMO (5′-ACCACACTAAGACATTGTAACAGAG). For Cox5aa experiments, a 5-bases mismatch MO was used as a control (5′-AGAGTCGAACGCCGGGTCGGAAGAT), otherwise the standard control MO of Genetools was used. Injected mRNAs were generated from linearized plasmids using the mMessage mMachine SP6 Kit (Ambion). To stimulate mitochondrial biogenesis, embryos were treated between 6–10 hours postfertilization (hpf) with 50 μM DCHC (dissolved in DMSO). Inhibition of respiration was achieved by treatment of embryos between 10–13 hpf with 0.1 μM rotenone (dissolved in DMSO, both from Sigma-Aldrich).

Cell culture. The following cell lines were used in this study: hTERT immortalized primary human fibroblast cell line 1BR3 (wild-type; University of Sussex), patient-derived cell lines MITO_HH0013 (carrying a homozygous pathogenic variant in MPV17) and MITO_HH009 (carrying a homozygous pathogenic variant in NFU1) (both collected at University Medical Center Hamburg-Eppendorf). Cells were cultured in MEMα containing 10% FBS and 1% penicillin/streptomycin (all from Life Technologies) at 5% CO₂ and 37°C. To obtain synchronous cultures in G0, cells were kept 2 to 3 days in starving medium (MEMα containing 0.1% FBS and 1% penicillin/streptomycin). Cells were treated with 0.5 μM Smoothened agonist (SAG; Enzo) during the last 24 hours of culturing to stimulate Hedgehog signaling activity. For immunofluorescence, cells were seeded on cover slips at a density of 1.7 × 10⁴ cells per cm², optionally synchronized in G0, and fixed with either 4% PFA or chilled methanol. KD of candidate genes was performed with custom made pools of 30 individual siRNAs all targeting an individual gene of interest or the negative control siPOOL 5 (all from siTOOLs Biotech): MTRR (NCBI gene ID 4552), ACOX1 (NCBI gene ID 51), OSBPL1A (NCBI gene ID 114876), CA5B (NCBI gene ID 11238), PHKA2 (NCBI gene ID 5256), TAZ (NCBI Gene ID 6901). Transfection of siRNA pools was carried out with Lipofectamine RNAiMAX (Thermo Fisher Scientific). Mitochondrial biogenesis was stimulated with 50 μM DCHC (dissolved in DMSO; Sigma-Aldrich). To this end, cells were treated with 50 μM DCHC for 48 hours after transfection; the last 24 hours were in starving medium if cells were to be synchronized in G0. Inhibition of respiration was achieved with rotenone treatment (100 nM). Treatment with MitoTempo (SML0737, Sigma-Aldrich) was carried out during synchronization with starving medium at a final concentration of 50 nM.

Mitochondrial DNA content. RNase A–treated genomic DNA from blood samples (heterotaxy patients and healthy subjects; Figure 1A) was diluted to a concentration of 5 ng/μl. Triplicate reactions were set up containing genomic DNA, iTaq Universal SYBR Green Supermix (Bio-rad), and one of the following primer pairs: (a) targeting the sequence of the mitochondrial-encoded 12s ribosomal RNA gene (Genebank ID 4549; F, 5′-CTAACCCCAGGGTTGGTCA; R, 5′-TACGCCGGCTTCTATTGACT) or (b) the exon-intron transition of exon 1 of the nuclear gene MLH1 (ENST00000231790.6; F, 5′-GTTCGTGGCAGGGGTTATT; R, 5′-TCCCTCCGTACCAGTTCTCA). Amplifications were carried out in an ABI PRISM 7300 Sequence Detection System (Applied Biosystems) or a LightCycler 480 Instrument II (Roche). Levels of the mitochondrial encoded 12s ribosomal RNA gene were compared with the nuclear encoded MLH1 gene to assess abundance of mitochondrial DNA.

ATP measurement. The ATPlite Luminescence ATP Detection Assay system (Perkin Elmer) was used for ATP measurements. Mitochondrial function or biogenesis was manipulated in zebrafish embryos and human fibroblasts as described above. Zebrafish embryos of the desired stage were collected and flash frozen in aliquots of 2 embryos. For measurements, embryos were placed on ice, immediately lysed in 30 μl mammalian cell lysis buffer (component of the kit), and transferred to a well of a black, Nunclon Delta Surface flat-bottom 96-well plate (Thermo Fisher Scientific). Human fibroblasts were plated in the same type of cell plate at a density of 1000 cells per well and were lysed by addition of 30 μl mammalian cell lysis buffer to 50 μl of medium and shaking at 400 rpm for 5 minutes at room temperature. To both types of lysates, 30 μl of substrate solution was added and shaken at 400 rpm for 5 minutes at room temperature. Wells containing no lysate or defined amounts of ATP were treated in parallel for blank controls as well as calculation of a standard curve. After 15 to 90 minutes of incubation in the dark, luminescence was measured using an infinite M200 plate reader and i-control 1.6 software (both from Tecan).

Measurement of oxygen consumption rate (OCR). The OCR was measured with a Seahorse XF96 device (Agilent). Zebrafish embryos were treated between 6–10 hpf with 50 μM DCHC (dissolved in DMSO) or vehicle. At tailbud stage, embryos were transferred to a XF96 plate (5 embryos per well). Basal oxygen consumption was measured 3 times over 2 minutes.

Whole mount in situ hybridization. Embryos were allowed to develop to the indicated stages. To prevent pigmentation, embryos older than 22 somites were treated with 0.003% 1-phenyl-2-thiourea (PTU, Sigma-Aldrich). At desired stages, the chorion was removed and embryos fixed in 4 % paraformaldehyde at 4°C overnight. In situ hybridizations were carried out following standard protocols (57). In situ probes were labelled with DIG by in vitro transcription from linearized plasmids encoding for the following mRNAs or parts thereof: charon (NM_212969), southpaw (spaw; NM_180967), cardiac myosin light chain 2 (cmlc2; AF425743.1), insulin (ins; NM_131056), 5-methyltetrahydrofolate-homocysteine methyltransferase reductase (mtrr; ENSDARG00000045398), mitochondrial ribosomal protein L38 (mrpl38; NM_212945.1), coenzyme Q9 homolog (coq9; NM_001098746.1), mitochondrial ribosomal protein S25 (mrps25; NM_001017590.1), coiled-coil-helix-coiled-coil-helix domain containing 3a (chchd3a; NM_214804.2), mitochondrial ribosomal protein L44 (mrpl44; NM_001030260.2), A kinase (PRKA) anchor protein 1b (akap1b; NM_001098179.2), acyl-CoA oxidase 1, palmitoyl (acox1; NM_001005933.2), oxysterol binding protein-like 1A (osbpl1a; BC162978.1), translocase of inner mitochondrial membrane 17 homolog B (yeast) (timm17b; NM_001113593.1), APEX nuclease (apurinic/apyrimidinic endonuclease) 2 (apex2; NM_200146.1), aminolevulinate, delta-, synthase 2 (alas2; NM_131682.2), solute carrier family 25 member 5 protein (slc25a5; AF506216.1), NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 1 (ndufa1; NM_001002153.3), isocitrate dehydrogenase 3 (NAD⁺) gamma (idh3g; NM_001243172.1), tafazzin (taz; NM_001001814.1), phosphorylase kinase, alpha 1b (phka1b; ENSDART00000160800.1), the closest homolog to human CA5B, carbonic anhydrase Va (ca5a; NM_001111201.1), aarF domain containing kinase 4 (adck4; ENSDART00000172195.1), and phosphorylase kinase, alpha 2 (phka2; ENSDART00000078202.4).

Cilia staining in cells and zebrafish embryos. Fixed G0-synchronized cells were permeabilized with PBS containing 0.1% Triton X-100 and then blocked with 3% BSA in PBS containing 0.1 % Tween 20 (PBST). To visualize cilia in the KV, fixed zebrafish embryos of specific stages were processed according to Jaffe et al. (58). Primary antibodies specific for acetylated tubulin (1:500 in embryos, 1:1000 on cells, clone 6-11b-1; Sigma-Aldrich), γTubulin (1:1000 in cells, T5192; Sigma-Aldrich) and COX IV (1:1000, clone 3E11; Cell Signaling Technology) were used. Fluorescent dye–labeled secondary antibodies were used to detect bound primary antibodies (Life Technologies). Stained embryos were dissected and flat-mounted between 2 cover slips in Vectashield medium containing Dapi (Vectorlabs). Cells were also mounted in Vectashield medium containing Dapi. Confocal stacks were generated with increments of 0.25 μm (cells) or 0.4 μm (zebrafish embryos) with sequential laser excitation and rendered into single-plane projections for analysis. The entire length of cilia was tracked with the help of an Intous3 6 × 8 USB tablet (Wacom) and measured with ImageJ.

qPCR. Total RNA was extracted from human fibroblasts or zebrafish embryos using the Quick-RNA MiniPrep kit (Zymo Research), which includes genomic DNA removal columns and on-column DNaseI digest of remaining genomic DNA. Reverse transcription was carried out using the ProtoScript II First Strand cDNA Synthesis Kit (New England Biolabs) or Superscript II Reverse Transcriptase (Thermo Fisher Scientific). Expression levels of specific mRNAs were assessed with the Universal Probe Library system (Roche), ABsolute QPCR Mix (with ROX, Thermo Fisher Scientific), and the following primers (and probes): TAZ (probe 68; F, 5′-TTCCCCTCATTCTCTGCTTG; R, 5′-CACTCCCACTTCTTCAGCTTG), PHKA2 (probe 25; F, 5′-TCGTGGATGGAGATGGTAGG; R, 5′-TTGGCGTTGAAAGAAGTGTG), OSBPL1A (probe 3; F, 5′-GAGTCTAAATCCAAGGTCATCTCG; R, 5′-GGACTGGAAAAGCTTTCAACC), MTRR (probe 2; F, 5′-TGGCCATATTCTTTGCATCTC; R, 5′-CAGCAGGTGGCGAGAATC), CA5B (probe 3; F, 5′-CAAATCTGACTGCGTTCCAA; R, 5′-ACACCGTGGACAGCAAATG), ACOX1 (probe 76; F, 5′-TCCTGGGTTTCAGGGTCATA; R, 5′-CCCAGACAGAGATGGGTCAT), SDHA (probe 80; F, 5′-GGACCTGGTTGTCTTTGGTC; R, 5′-CCAGCGTTTGGTTTAATTGG), GLI1 (probe 7; F, 5′-CCAGCCAGAGAGACCAACAG; R, 5′-CCCGCTTCTTGGTCAACTT), acox1 (probe 3; F, 5′-TCTGAACATCGCCCAGGTG; R, 5′-GGCATATGCAGTAGCCAGAAG), caVa (probe 51; F, 5′-TCCGCAAACTATCAATCCTCA; R, 5′-TCAATAGGTGACTGTCGATCTCC), mtrr (probe 17; F, 5′-TGCCGTCACAAAGACAAAGA; R, 5′-CTCAGGGTCCCATTATGCAC), osbpl1a (probe 82; F, 5′-CTGTTGGAGGACCGAGTCA; R, 5′-CTGGCATCCCTGAGGAAC), phka2 (probe 159; F, 5′-GGCACTTACACTTTTACTCCACAG; R, 5′-TGGTTATCCAGAGCCAAGTAGAA), taz (probe 62; F, 5′-GGCAGTTGTGGAACCTTAACA; R, 5′-GGAGTGAAATTCTCTTGTGAAACA), gapdh (probe 147; F, 5′-CAGGCATAATGGTTAAAGTTGGTA; R, 5′-CATGTAATCAAGGTCAATGAATGG). Reactions were run in triplicate in a LightCycler 480 Instrument II and results were analyzed with LightCycler 480 Software, version1.5 (both from Roche).

Exome sequencing and variant analyses. Whole-exome sequencing (WES) was performed on an Illumina HiSeq2500 sequencer using the Nimblegen SeqCap EZ Human Exome kit at the CCHMC DNA Sequencing and Genotyping core. Data are deposited at the NCBI Database of Genotypes and Phenotypes (dbGAP) under the reference phs001814.v1.p1. Alignment, preprocessing, and variant calling was performed using the Genome Analysis Toolkit (GATK) at the CCHMC Bioinformatics Core. Variant filtering and annotation was done using Golden Helix SVS software. For the initial variant assessment and definition of mitochondrial disease incidence, genes were annotated against RefSeq 105 v2 and restricted to homozygous and X-linked variants present in exon regions, which resulted in protein alterations including missense, frameshift, splice site, stop loss, stop gain, and small in-frame insertions and deletions. Only variants with a read depth of at least 10 were analyzed. A minor allele frequency (MAF) of less than 1% was used to further filter rare variants. Population allele frequencies were obtained from the Database of Single Nucleotide Polymorphisms (dbSNP) build 147, the National Heart, Lung, and Blood Institute (NHLBI; ESP6500SI-V2-SSA137), and the Exome Aggregation Consortium (ExAC). A high priority list of 1174 genes that localize to mitochondria or cause mitochondrial disorders was investigated (according to the Human MitoCarta2.0 [ref. 36] and the BCM-MitomeNGS [Baylor Genetics]). These were further prioritized by Combined Annotation Dependent Depletion (CADD) and dbNSFP Functional Prediction programs SIFT, Polyphen2 HVAR, MutationTaster, MutationAssessor, FATHMM, and FATHMM MKL Coding. For the variant burden analyses, we applied additional quality filters to the heterotaxy cohort and controls from the 1000 Genomes study. Variants in the heterotaxy cohort that passed a read depth of at least 10, genotype quality of at least 30, call rate of at least 0.8, and a cohort allele frequency no greater than 0.2 were used. For the 1000 Genomes participants, the 1000 Genomes Pilot mask obtained from The International Genome Sample Resource (IGSR) was applied. Variants in both cohorts were further restricted to regions covered by the Nimblegen Seq capture platform. Missense variants with a MAF below 0.1% in ExAC were filtered, followed by MetaSVM analysis to predict damaging and tolerated variants. We normalized damaging variants to tolerated variants to obtain the overall damaging variant burden ratio in mitochondria-associated genes and control gene lists. Control gene lists were created by identifying genes from external databases and clinical panels, including the Candidate Cancer Gene Database for cancer genes, Autism/ID Xpanded Panel from GeneDx for autism genes, and Orphanet, Online Mendelian Inheritance in Man (OMIM) Morbid, and The Developmental Disorders Genotype-Phenotype Database (DDG2P) for random disease genes (Supplemental Table 4). Control gene lists were matched to the mitochondria-associated gene list by total exon length (4,969,960 base pairs in mitochondria genes) with differences between the gene lists less than 1000 base pairs.

Verification of MO efficiency. MO efficiency was either verified by antibody staining (Phka2: rabbit anti-PHKA2, Proteintech, catalog no. 24658-1-AP, lot 20988, 1:100; Osbpl1a: rabbit-anti-Osbpl1a, GeneTex, catalog no. GTX107722, clone C3, 1:200) or by RT-PCR to show splice blocking. Injection of the Mtrr splMO resulted in deletion of the start codon as verified with the following primers: 5′-GATTACCTGTTCCTCCCGTTC and 5′-ATGAGAAAGATGGTCGTGTGG. The Taz splMO interfered with regular splicing of exon 2, which was tested with the following primers: 5′-CGGAAGTGTGGTTTAATGCTG and 5′-TGATAAGGGGTGTGTCCTGAG. See complete unedited gels in the supplemental material.

Cloning and in vitro synthesis of RNA. Murine Ppargc1a (Pgc1a) (NM_008904.) was amplified from mouse cDNA generated from total RNA isolated out of lineage-negative hematopoietic progenitor cells and cloned via BamH/XhoI into pCS2⁺. A clone of human MTRR was purchased from Dharmacon (clone 40007734) and was subcloned by directional TOPO cloning into pcDNA3.1 before positional mutagenesis to induce the R552Q mutation as detected in heterotaxy patients. Similarly, the open reading frame of human TAZ (NM_000116) was amplified from cDNA of human fibroblasts, cloned via ClaI/StuI into pCS2⁺, and the R253Q mutation was introduced. For injection into zebrafish eggs, plasmids were linearized using NotI (pCS2⁺-mPgc1a) or PmeI (pcDNA3.1 constructs), respectively. Capped RNA was transcribed using either the SP6 or T7 mMessage mMachine Kit (Ambion) or the AmpliCap SP6 High Yield Message Maker Kit (Cellscript). In preparation for in situ hybridization experiments, fragments of open reading frames of candidate genes identified in the heterotaxy exomes (see also “Whole-mount in situ hybridization”) were amplified using either GoTaq (Promega) or the Expand Long Range polymerase (Roche) and cloned into pCRII by TOPO TA cloning before linearization and in vitro transcription of DIG-labeled antisense RNA using SP6 or T7 RNA polymerases (New England Biolabs).

Imaging. Confocal images were acquired on a TCS SP5II with LAS AF software (both Leica). Cells were analyzed using an AxioPhot (Zeiss), an ORCA-03G camera, and HoKaWo software (both Hamamatsu). Whole-mount embryos were imaged using a Leica M125 equipped with a Leica IC80 HD camera. Movement trajectories of Tetrahymena thermophila were analyzed with a Discovery V8 equipped with an Axiocam 506 mono and ZEN software (all Zeiss).

Measurement of ciliary beat frequency (CBF) and KV flow. Embryos were mounted in 1% agarose and set under the ×100/1.30 NA oil immersion objective lens on a Nikon Eclipse Ti-U inverted microscope at room temperature (24°C–25°C). All images were taken with the dorsal roof of the KV facing the objective lens. Bright field images were recorded with a FASTCAM MC2 camera (Photron Europe Limited) controlled with Photron FASTCAM Viewer software. Native particles were filmed at 60 frames per second (fps) for 30 seconds whereas cilia were recorded at 500 fps for 2 seconds. KV flow and CBF measurements were analyzed using Fiji software as previously described (23).

Analysis of Tetrahymena thermophila. The wild-type strain CU428 (gift from Dorota Wloga, Nencki Institute, Warsaw, Poland) was maintained in SPP medium (1% proteose-peptone, 0.2% D⁺-glucose, 0.1% yeast extract, 0.03 g/l EDTA-Fe³⁺ salt; all from Sigma-Aldrich) in the dark and passaged every month. Exponentially growing cells at a density of 2 × 10⁵ to 4 × 10⁵ cells per ml were treated for 1 hour at 30°C with 40 rpm shaking with 50 μM DCHC, 180 μM rotenone, or DMSO as control, to manipulate mitochondrial function. To test phagocytosis efficiency, 100 μl cell suspension was supplemented with 10 μl drawing ink 17 (Pelikan) and incubated without shaking for 10 minutes at 30°C. Cells were then fixed with 4% PFA and vacuoles of at least 100 cells from 5–6 independent experiments were counted using an Axiophot microscope (Zeiss). To measure swimming behavior, 1.2 ml of pretreated cells were transferred to a new 6-cm dish. In equal darkfield settings for all treatments, cells were imaged with an exposure time of 10 seconds. Trajectories of over 400 cells from 2 independent experiments were measured for each treatment with ImageJ.

Cryo-electron microscopy. Scanning transmission electron microscopy (STEM) tomography was performed as previously described with some modifications (59). Sections (700-nm thick) were cut from the epoxy resin block parallel to the plane of the sapphire disc with an Ultracut UCT ultramicrotome (Leica) equipped with a diamond knife (Diatome). Before mounting the slice onto a special copper grid for tomography with parallel bars, the grid was plasma-cleaned with an Edwards plasma cleaning system. A droplet of 10% (wt/vol) poly-L-lysine (Sigma-Aldrich) in water was added on the sample and dried for 5 minutes at 37°C. The grid with the slice on it was treated on both sides with 15 μl of a solution containing 25 nm gold particles (Aurion) diluted 1:1 with water. Finally, the grid was coated on both sides with a 5 nm carbon layer using a BAF 300 electron beam evaporation device (Balzers). Images were recorded in the STEM mode with a Jeol JEM-2100F (Jeol) with an acceleration voltage of 200 kV and a bright field detector (Jeol) at a size of 1024 × 1024 pixels and an exposure time of 22 seconds. Images were recorded at different tilt angles ranging from −70° to +70° with a tilt increment of 1.5°. The complete tilt series contained 97 individual images. The 3D stack was reconstructed by a weighted back-projection algorithm with the help of the IMOD software (60). For segmentation the software Amira (FEI) was used.

Statistics. Statistical analysis was carried out with Graphpad Prism 6, 7 and JMP Pro 14. N numbers, P values, and applied tests are indicated in figure legends. P values less than 0.05 were considered significant. Samples were tested for Gaussian distribution and equality of SD. Generally, 2-tailed tests were applied. However, 1-tailed tests were applied as indicated in figure legends, when biological changes were expected to be unidirectional, such as correct versus wrong heart looping.

Study approval. The study was performed in accordance with the Declaration of Helsinki protocols. Prior to inclusion in the study patients or their respective parents gave their written informed consent. The local ethics committees in Berlin, Ulm, and Hamburg as well as the CCHMC Institutional Review Board approved sample collection, distribution, and analysis. The Veterinary Care Unit at Ulm University and the animal welfare commissioner of the regional board for scientific animal experiments in Tübingen, Germany, approved husbandry as well as manipulation of zebrafish. Zebrafish experiments were performed according to the European Union Directive 86/609/EEC for the protection of animals used for experimental and other scientific purposes.