Rare coding variants in axonemal dynein heavy chain genes are associated with adolescent idiopathic scoliosis

Adolescent idiopathic scoliosis (AIS) is the most common pediatric musculoskeletal disorder worldwide, characterized by atypical spine curvatures in otherwise healthy children. Human genetic studies have identied candidate genes associated with AIS, however, only a few of these genes have been shown to recapitulate adult-viable scoliosis in animal models. To further dene susceptibility loci for AIS, we performed whole exome sequencing on a cohort of 195 Han Chinese AIS patients and 229 healthy controls. We identied members of the axonemal dynein family associated with both sporadic and familial AIS. We demonstrate that disruption of the dynein axonemal heavy chain 10 (dnah10) gene results in recessive adult-viable scoliosis in zebrash. These dnah10 mutant zebrash display reduced ependymal cilia beating and a disassembly of the Reissner ber in the hindbrain and spinal canal, concurrent with the onset of body curvatures. Altogther, these results demonstrate that mutations in axonemal dynein genes are linked with human AIS and suggest that ependymal cell cilia function plays an essential role in maintaining spine alignment in humans. Using CRISPR/Cas9 targeted mutagenesis, we functionally veried that a frameshift-truncating mutation in the dynein axonemal heavy chain 10 (dnah10) gene resulted in adult-viable scoliosis in zebrash. High-speed time-lapse confocal imaging of motile cilia within the spinal cord central canal demonstrated alterations in motile cilia beating. Using an endogenous RF-labeling scospondin-GFP ut24 knock-in zebrash line, we demonstrated a range of RF assembly defects concurrent with the onset and severity of spine curvatures in dnah10 ut28 mutant zebrash. To our knowledge, this is the rst study to associate rare variants of several axonemal dynein gene with human AIS, with relevant functional testing and demonstration that dnah10 is essential component of motile cilia dynamics and spine morphogenesis in zebrash. This study provides key evidence that alterations in ependymal cell motile cilia function may underlie AIS in humans.


Introduction
Adolescent idiopathic scoliosis (AIS) is a common pediatric disorder characterized by atypical three-dimensional curvatures, with rotation of the spinal column (Cheng, Castelein et al. 2015). The prevalence of AIS is approximately 1.0-5.0% among children worldwide (Zhang, Guo et al. 2015, Yılmaz, Zateri et al. 2020). Approximately 10% of AIS patients suffer severe physical and mental issues due to progressive spine curvature (Cheng, Castelein et al. 2015). Analysis of familial scoliosis and the recognition of a high concomitance of scoliosis in twin-studies demonstrate that AIS is a multifactorial disease with a strong genetic predisposition (Hadley Miller 2000, Wise, Gao et al. 2008. Hypothesis-free genome wide association study (GWAS) analyses from different ethnic backgrounds have mapped a range of candidate loci and implicated several candidate genes associated with this disorder, including LBX1, BNC2, GPR126, AJAP1 and PAX1 (Takahashi, Kou et al. 2011, Kou, Takahashi et al. 2013, Ogura, Kou et al. 2015, Sharma, Londono et al. 2015, Zhu, Xu et al. 2017). Recent, exome sequencing studies has identi ed additional rare variants associated with AIS in different genes, including within POC5, HSPG2, AKAP2, SLC39A8 and FLNB (Baschal, Wethey et al. 2014, Patten, Margaritte-Jeannin et al. 2015, Li, Li et al. 2016, Haller, McCall et al. 2018, Jiang, Liang et al. 2020. Additional analysis of exome data have identi ed an enrichment of potentially deleterious variants in FBN1/2, extracellular matrix (ECM) and cilia-associated genes in AIS patients (Buchan, Alvarado et al. 2014, Haller, Alvarado et al. 2016, Baschal, Terhune et al. 2018. Overall, these data suggest that AIS is a complex polygenic disease with signi cant genetic heterogeneity. Therefore, additional genetic studies are needed to expand our understanding of the genomic architecture of AIS. alterations in cerebrospinal uid (CSF) ow, disruption of the Reissner ber (RF), and urotensin neuropeptide signaling to the dorsal musculature are antecedents of scoliosis in zebra sh (Bagnat and Gray 2020).
The RF, is a glycoprotein-rich extracellular thread, which stretches from the third cerebral ventricle down the central canal in most chordate lineages (Munoz, Kahne et al. 2019). Dysfunction of motile cilia in the central canal of the spine leads to disrupted RF assembly and severe ventral ward body curvatures (Cantaut-Belarif, Sternberg et al. 2018). Importantly, the continuous assembly of the RF during larval stages of development in zebra sh is essential for maintaining spine stability in larval zebra sh (Rose, Pompili et al. 2020, Troutwine, Gontarz et al. 2020) and the assembly of the RF appears to require the function of motile cilia components, including kinesin family member 6 (kif6) (Troutwine, Gontarz et al. 2020) and protein kinase 7a (ptk7a) (Rose, Pompili et al. 2020). Finally, loss of CSF ow and RF assembly can alter the expression of urotensin II-related neuropeptides in the central canal (Zhang, Jia et al. 2018, Vesque, Anselme et al. 2019, Lu, Shagirova et al. 2020, which has been shown to signal to the dorsal musculature to help maintain a straight body axis in larval zebra sh (Zhang, Jia et al. 2018). Altogether, analysis of scoliosis in zebra sh have begun to elucidate a robust pathway for control of spine morphogenesis. However, there remains a lack of evidence linking the pathogenesis of scoliosis in these zebra sh models to genetic studies of human AIS.
In this study, we report a whole exome sequencing (WES) on 195 sporadic AIS patients and 229 in-house controls to identify novel AIS-associated mutations enriched with variants predicted to be deleterious. We inferred variants in ten genes encoding a variety of dynein axonemal heavy chain proteins. Using CRISPR/Cas9 targeted mutagenesis, we functionally veri ed that a frameshift-truncating mutation in the dynein axonemal heavy chain 10 (dnah10) gene resulted in adult-viable scoliosis in zebra sh. High-speed time-lapse confocal imaging of motile cilia within the spinal cord central canal demonstrated alterations in motile cilia beating. Using an endogenous RF-labeling scospondin-GFP ut24 knock-in zebra sh line, we demonstrated a range of RF assembly defects concurrent with the onset and severity of spine curvatures in dnah10 ut28 mutant zebra sh. To our knowledge, this is the rst study to associate rare variants of several axonemal dynein gene with human AIS, with relevant functional testing and demonstration that dnah10 is essential component of motile cilia dynamics and spine morphogenesis in zebra sh. This study provides key evidence that alterations in ependymal cell motile cilia function may underlie AIS in humans.

AIS and control subjects
This study was approved by the Ethics Committee of Xiangya Hospital of Central South University (Changsha, China) and all participants and their legal guardians involved in this project have signed a written consent form. A total of 195 unrelated Han Chinese AIS patients were recruited from Xiangya Hospital. The in-house control cohort consisted of 229 healthy Han Chinese subjects from The Han Chinese Genomes Database (PGG.Han), a population genome database to serve as the central repository for the genomic data of the Han Chinese genome database (https://www.hanchinesegenomes.org/HCGD/index) (Gao, Zhang et al. 2020). The inclusion criteria for AIS were as follows: (1) AIS diagnosis for each case was con rmed by three certi ed senior spine surgeons.
Patients with congenital vertebral malformations, connective tissue diseases, skeletal developmental delay or other scoliosis related syndromes were excluded; (2) Cobb angle > 10°; (3) No second-degree family members of AIS cases were diagnosed as scoliosis or inherited diseases. The mean Cobb angle of these AIS patients was 41.88 ± 18.35°. The female-male ratio was 4 (156 females vs. 39 males), which following the natural sex distribution of AIS population (Table S1).
The second exome sequencing cohort consists of 28 IS families enrolled as previously described (Baschal, Wethey et al. 2014). A diagnosis of IS required a standing anteroposterior spinal radiograph showing ≥ 10°c urvature by the Cobb method with pedicle rotation, and no congenital deformity or other co-existing musculoskeletal disorder. Pedigrees were created using PedigreeXP.
Whole exome sequencing and data preprocessing Genomic DNA was extracted from the venous peripheral blood samples of AIS cases. Exome sequencing within sporadic AIS patients was performed on BGISEQ-500 sequencing platforms (BGI Group, Wuhan, China) according to the manufacturer's protocol. After high-throughput sequencing, the raw data of each case were generated and stored in the FASTQ format. Clean data were kept by removing reads including sequencing adapters, low-quality (base quality ≤ 20) and uncertain base (N > 10%). We mapped clean reads to the human reference genome (hg19/GRCh37) using the Burrows-Wheeler Aligner (BWA; Version 0.7.10) and removed the duplicate reads marked by Picard tools (Version 1.117). Then we used the Genome Analysis Toolkit (GATK; Version 3.6) to perform local realignment, Base Quality Score Recalibration (BQSR) and variants calling (Using HaplotypeCaller and output in a genomic variant call format [GVCF] le).
Exome sequencing within 28 IS families was conducted and ltered as described previously. Brie y, 3-6 affected individuals were sequenced per family using the Illumina TruSeq exome kit and sequenced with a 2 x 100 bp run on the Illumina HiSeq 2000 at the University of Colorado Denver Genomics and Microarray Core facility. Reads were aligned to the Hg38 genome and variants were identi ed with FreeBayes.

Joint calling, quality control and variants ltration
We performed joint calling in all sporadic samples (195 AIS and 229 controls) using GATK. To generate an accurate calling result, We performed variant-level quality control (QC). we ltered out the variants with low coverage due to uneven coverage of WES by remove variants which missing rate more than 20% in AIS. To further lter out low-con dence variants, Variant Quality Score Recalibration (VQSR) from GATK was performed. In this step, we select the ltering indicators including read depth, strand bias, mapping quallity, read position bias.To further identify the AIS-related candidate variants, we used the following lters: (1) exclude high frequency SNPs (MAF ≥ 0.01) in an in-house control dataset of Han Chinese, Han Chinese in Beijing (CHB) and Southern Han Chinese (CHS), from 1000 Genomes Project (1000G) dataset; (2) Variants with highest frequency in AIS cohort among AIS, in-house control, CHB and CHS cohorts were reserved; and (3) Variants classi ed as "HIGH" or "MODERATE" impact by Variant Effect Predictor (VEP) were kept. All the ltered deleterious candidate variants were subjected to gene-burden analysis.
For the exome data from 28 IS families. Variants were retained that were present in all sequenced family members, had an allele frequency of £ 5% within ExAc, and were scored as damaging by at least one predictive functional algorithm (SIFT, Polyphen2, LRT, MutationTaser).
Sanger sequencing DNA from additional affected and unaffected family members was Sanger sequenced across the variants of interest detected through exome sequencing. Primers were designed using Primer3 and obtained from Integrated DNA Technologies. PCR was conducted in 20 µL reactions with 10 µL Premix D (Epicentre Biotechnologies), 0.2 µL Taq Polymerase (Sigma), 20 ng genomic DNA, and 20 µM Forward and Reverse Primers. PCR reactions were run with a touchdown PCR protocol on a SimpliAmp Thermocycler (Fisher Scienti c). Annealing temperatures were set with the delta function to begin at 65°C and decrease by 0.3 °C each cycle, ending at 55 °C. PCR reactions were cleaned up for Sanger sequencing using ExoSAP-IT Express (Thermo Fisher). Sanger sequencing was performed by Quintara Biosciences and chromatograms were analyzed using the CodonCode Aligner v9.0 (CodonCode Corporation, https://www.codoncode.com/index.htm).

Gene-burden analysis and functional analyses
To perform gene-burden analysis, we determined two indicators in gene-level, The samples carrying proportion of harmful variation and deleterious variant count per sample, in AIS cases and in-house controls by the candidate variants. Considering sequence discrepancy may lead to loci missing, we corrected the sample size for each gene using the equation below.
Where m denotes number of variants, n indicates observable sample size, and AF represents derived allele frequency. The gene-level samples carrying proportion of harmful variation generated by dividing the corrected sample size by counting samples that contain deleterious variants for each gene. Gene-level deleterious variant count per sample was calculated as the number of deleterious variants in one gene divide corrected sample size in this gene. To avoid the bias of gene length, we correct these two indicators by dividing them by length. A twosided Fisher's exact test was performed with deleterious variant count per sample in two datasets. The difference was estimated between deleterious variant proportions in the two datasets. Further, we used a systematic functional annotation approach, including Gene Ontology (GO) analysis, KEGG pathway, gene family enrichment analysis and Protein Annotation Through Evolutionary Relationship (PANTHER) classi cation system to interpret these AIS-associated genes (Mi, Muruganujan et al. 2019). Gene family annotation information was based on HGNC database (HUGO Gene Nomenclature Committee).

Zebra sh husbandry
All in vivo experiments were performed following approved Institutional Animal Care and Use Committee (IACUC) protocols. All zebra sh procedures were approved by the Animal Studies Committee at the University of Texas at Austin (AUP-2018-00342). Animals were maintained and raised at 28.5 °C with 14/10 hours of light-dark periodicity.

Skeletal staining
Animals were xed in 10% formalin overnight, then incubated in 100% acetone overnight, and stained in Bone/Cartilage Stain (59.5% Ethanol, 5% Acetic Acid, 0.015% Alcian Blue and 0.005% Alizarin Red) solution at 37°C overnight. Solution and soft tissue were cleared with 1% KOH until only stained skeletal tissue was left. After clearing, specimens were further moved to 25%, 50% and 80% glycerol. Zebra sh and stained skeletal specimen were imaged using an Olympus SZX2-ILLT microscope (Olympus, Tokyo, Japan).

CRISPR-mediated zebra sh F0 screen
To rapidly test whether these AIS candidate genes may be involved in spine development, we used a CRISPR/CAS9 screening approach by injecting four guide RNAs to redundantly target the same genetic loci to generate thorough gene disruption (Wu, Lam et al. 2018). In brief, sgRNAs were synthesized in vitro with T7 RNA Polymerase (New England Biolabs, Ipswich, MA, USA) and puri ed by Zymo RNA Clean & Concentrator kit (Zymo, Irvine, CA, USA). A mixture of 5μM of Cas9 protein (Integrated DNA Technologies, Coralville, IA, USA) and 31 μM of sgRNA (25 μM each) was made and incubated for 5 minutes at room temperature. About 500 pL of the mixture of Cas9 and sgRNAs was injected into the yolk of 1-cell stage zebra sh embryos. Embryos were raised and screened at 3-and 5-days post fertilization (dpf). A total of six genes of interest, including dnah1, dnah3, dnah7l, dnah7, dnah9 and dnah10, were screened using this method. Oligos and gRNA target sequences are listed in Table S2.

Identi cation of dnah10 mutation
To preserve and isolate the dnah10 mutation, we outcrossed the F0 scoliosis zebra sh with AB wild-type sh to generate F1 zebra sh. Then we in-crossed the F1 zebra sh and generated scoliosis F2 zebra sh. One scoliosis F2 zebra sh was chosen and outcrossed with wild-type and other transgenic zebra sh for further mutation identi cation and functional analyses. To determine the mutation of the scoliosis dnah10 mutant, genomic DNA from heterozygous F3 (scoliosis F2 outcross with wildtype) zebra sh was extracted from tissue samples using 50 mM NaOH according to the protocol. DNA segments around four gRNA targets in the dnah10 gene were ampli ed and Sanger sequenced (PCR primers listed in Table S3). Variants were identi ed using ApE (Version 0.55) and annotated using SnpEff (Version 4.3T).

Whole mount in situ hybridization
Whole mount in situ hybridization was performed to detect dnah10 expression according to a previously published protocol (Thisse and Thisse 2008). For the generation of dnah10 probes, a ~1.2kb fragment of the dnah10 transcript (ENSDART00000087708.6) was ampli ed from cDNA of AB wild-type (PCR primers see Table  S3). The fragment of dnah10 was then cloned into a pGEM-T Easy Vector (Promega) and validated by Sanger sequencing. The template of dnah10 was generated by ApaI restriction enzyme cutting of the recombinant plasmid. The antisense DIG-labeled RNA probe of dnah10 was synthesized with T7 RNA polymerase.

In vivo imaging and immuno uorescence
The transgenic zebra sh strains scospondin-GFP ut24 and Tg (Foxj1:Arl13b-GFP) were created as previously described (Konjikusic, Yeetong et al. 2018, Troutwine, Gontarz et al. 2020. Tg(Foxj1:Arl13b-GFP); dnah10 ut28 and scospondin-GFP ut24 ;dnah10 ut28 lines were generated by cross homozygous dnah10 ut28 to each transgenic zebra sh and incrossing of their progeny to generate the homozygote. For live imaging, larvae were anaesthetized in 0.16% tricaine and embedded in 1% low-melt agarose. Time-lapse videos (from 10 s to 10 mins) were acquired at 0.1 to 33 Hz imaging frequency. Videos were cropped and edited by Adobe Premiere Pro (Version 12.0). Reissner ber and oor plate immuno uorescence staining was performed as previously described (Troutwine, Gontarz et al. 2020). In brief, samples were xed, permeabilized and incubated with primary antibodies. After washing and incubation with secondary antibodies, specimens were xed and lmed immediately. Confocal images and videos were taken using a Nikon Ti2E / CSU-W1 spinning disc confocal system. Images were analyzed and reconstructed using Fiji (Schindelin, Arganda-Carreras et al. 2012). ROIs were manually adjusted on max intensity Z-projection.

Ciliary movement analysis
All the videos for ciliary movement analysis were acquired at 33 Hz frequency. The range of ciliary motion was observed in a maximal intensity Z projection of 10 s time-lapse confocal video. Cilia beat frequency was calculated by Fiji (Schindelin, Arganda-Carreras et al. 2012). We drew a one-pixel width straight line across the cilia to create kymographs to obtain the time trajectory of cilia ( Figure S5). Kymographs generated a wave diagram for each cilium, and the beat frequency was calculated by pixels for one second (33 pixels for 33 Hz video) divided by pixels for one beat cycle. Cilia were sorted into "immobile cilia" and "motile cilia" categories according to their kymograph result ( Figure S5).

Statistical analysis
Statistical methods used for the gene-burden analysis was described previously. The ciliary data analyses were performed using Fiji (ImageJ) and GraphPad Prism (Version 8.0.2). Changes between the groups were examined using Student's t test. P < 0.05 was considered signi cant.

Results
Identi cation of loci associated with AIS in humans using gene-burden analysis To de ne the biological function of our candidate AIS susceptibility loci, we performed enrichment analysis on our set of 256 inferred genes. The most enriched gene family was "dynein, axonemal" (P = 3.66 × 10 − 10 ) ( Fig. 1A; Table 1, S5), implicating 10 genes associated with AIS in our cohort. An orthogonal enrichment analysis showed that genes encoding microtubule-binding motor proteins, were the most enriched protein class in the 256 candidate genes ( Fig. 1B; Table S6). Finally, KEGG pathway enrichment analysis of our candidate genes showed the most signi cance for the "ECM-receptor interaction" (P = 7.49 × 10 − 14 ) ( Fig. 1C; Table S7). Intersection of all three gene enrichment analyses showed an overlap of 10 dynein axonemal heavy chain genes common to Protein Class Ontogeny and Gene family enrichment analyses (Fig. 1D); however, there were no common candidate genes identi ed in the KEGG pathway enrichment analysis (Fig. 1D). We next analyzed DNAH variants detected in a second exome sequencing cohort of 28 IS families, ve families of which were previously published (Baschal, Terhune et al. 2018). Low frequency variants (minor allele frequency ≤ 5%) passing ltering were present within all sequenced members of the family (3-6 individuals sequenced per family) and were predicted to have damaging functional consequences (see Methods). Two variants within DNAH1 and DNAH6 were detected within two families, and one variant in DNAH7, DNAH8, DNAH14 and DNAH17 were each detected within one family. All variants were missense except for DNAH7, which was a splice site variant. One family possessed variants in both DNAH1 and DNAH6 (Family A) (Table S9).
We next performed Sanger sequencing of each variant within additional family members affected and unaffected with IS for which DNA was available. Pedigrees and corresponding genotypes are provided in Fig. 2. Of 8 variants, 6 variants segregated with the IS phenotype within the families. The variants that did not validate were within DNAH7 (Family E) and DNAH17 (Family F). The remaining variants segregated with the scoliosis phenotype, with the exception of an unaffected individual in Family A (II:9) who possessed the DNAH6 variant of interest and an affected individual who did not possess the DNAH8 variant of interest (III:6). Unaffected males III:1 in Family B and III:1 in Family D are obligate carriers who possessed the variants of interest, but are expected to as they have children affected with scoliosis.

A C-terminal truncating mutation of dnah10 causes scoliosis in zebra sh
Because several dynein axonemal gene family were enriched in our cohort of human AIS patients, we next sought to functionally test which axonemal dynein genes were involved in spine stability in zebra sh (Fig. 1D). Of these, six zebra sh dnah genes, with clear human homologues, were screened using a multiplexing CRISPR/Cas9 approach (Wu, Lam et al. 2018), employing four independent gene speci c guide RNAs injected with Cas9 protein (Table S2). Systematic screen of phenotypes at 3-, 5-days post fertilization (dpf) showed that targeting mutagenesis to dnah1 and dnah10 resulted in viable larval sh with atypical body curvatures (Fig. 3A, B, Table  S10). However, at 30 dpf only dnah10 targeted F0 zebra sh displayed adult-viable scoliosis (Fig. 3C).
We next sought to generate a stable genetic mutant of dnah10 by targeted CRISPR/Cas9 mutagenesis to the region encoding the microtubule binding domain of the Dnah10 protein. To identify a stable mutant line, we outcrossed F0 CRISPR/Cas9-injected zebra sh with scoliosis to expediently isolate deleterious dnah10 mutations in F1 families. We next incrossed siblings within distinct F1 families to screen for heterozygous carriers as shown in the schematic (Fig. 3D). Progeny from one of these crosses showed recessive scoliosis in adult zebra sh (Fig. 3E). Finally, we bottlenecked a single mutation from an individual F2 male dnah10 mutant zebra sh. From this progeny we identi ed a single homozygous frameshift variant in dnah10 (c.10304_10307del [p. N3435SfsX20]) that is predicted to generate a frameshift and premature stop codon 20-amino acids after N3435 in exon 60 of the dnah10 loci (Fig. 3F, S3), which was assigned a speci c allele dnah10 ut28 . Longitudinal analysis of the progression of body axis curvature phenotypes from crossing homozygous dnah10 ut28/ut28 mutants showed the penetrance of scoliosis was approximately 79% at 30 days post fertilization (dpf) (n = 165; Fig. 3G).
To determine the dnah10 expression pattern at different developmental stages, we cloned a ~ 1.2 kb piece of the dnah10 (ENSDART00000087708.6) cDNA to generate an immunolabeled riboprobe and assayed at a variety of developmental stages. Whole-mount in situ hybridization showed dnah10 expression at the 8-10 somite stage (12 hours post fertilization [hpf]) in a ubiquitous expression pattern. At 24 hpf, dnah10 expression was diffuse throughout the body but with more robust staining within the brain and pronephric tubules. After 48 and 72 hpf, dnah10 expression was found to be much more localized within the brain and pronephric tubules (Fig. 3H).
Dnah10 function is required for ependymal cilia motility in the ventricular system Dnah10 (called DHC1 in Chlamydomonas) functions as an inner arm dynein motor of motile cilia, leading to reduced motility when mutated in Chlamydomonas (Myster, Knott et al. 1997, Kamiya andYagi 2014). DNAH10 knockdown in the Trypanosoma brucei agellate caused severe motility defects of cilia, without obvious ultrastructural defects of the axoneme (Zukas, Chang et al. 2012). To determine whether dnah10 ut28/ut28 mutant zebra sh show alterations in the development or motility of motile cilia in the central canal we crossed it to a transgenic Tg(Foxj1a:Arl13B-GFP) ut22 strain, which labels motile ependymal cilia in the brain and central canal of zebra sh (Konjikusic, Yeetong et al. 2018), without affecting the body axis (Fig. 4A). Analysis of the central canal in fresh-xed animals showed no difference in the density of Arl13B-GFP labeled foxj1a expressing motile cilia comparing wild-type (WT) (Fig. 4C, E) and homozygous dnah10 ut28/ut28 mutant transgenic sh (Fig. 4D, E) (67.1 ± 5.0 vs 66.5 ± 5.3 cilia in 200 µm, n = 15).
Confocal imaging using slow acquisition speed of Arl13B-labeled motile cilia in the central canal of living zebra sh produce a characteristic fan-like appearance, as the result of blurring demonstrating the range and direction of ciliary motion (Borovina, Superina et al. 2010). This same fan-like shape was observed in WT Tg(Foxj1a:Arl13B-GFP) ut22 zebra sh (Fig. 4G). In contrast, dnah10 ut28/ut28 ; Tg(Foxj1a:Arl13B-GFP) ut22/+ transgenic mutant zebra sh showed fewer fan-shaped cilia (Fig. 4H), which establishes that cilia motility in dnah10 ut28/ut28 mutants is reduced. Next, we analyzed of cilia motility using high-speed (33 Hz) spinning disc confocal time-lapse imaging, which consistently revealed rapid beating of cilia in WT Tg(Foxj1a:Arl13b-GFP) transgenic animals (Video S1). In contrast, dnah10 ut28/ut28 ; Tg(Foxj1a: Arl13B -GFP) ut22/+ transgenic mutant zebra sh displayed a consistent reduction in central canal cilia motility (Video S2). Kymographic time trajectory analysis of several WT Tg(Foxj1a:Arl13b-GFP) motile cilia showed a distinctive waveform pattern for central canal cilia, while dnah10 ut28/ut28 ; Tg(Foxj1a: Arl13B -GFP) ut22/+ mutants displayed severely reduced waveform ( Fig. 4I). Interestingly, we observed that the loss of central canal cilia motility in dnah10 ut28/ut28 mutants was at times regionally variable within the same animal at 3 dpf. For example, we would sometimes observe reduced or absent central canal cilia motility in a more anterior region, while the more posterior region could display typical motility, in the same individuals (Fig. 4J, J' and Video S3). To assess which region of the ventricular system was most consistently susceptible to a reduction of motile cilia defects, we imaged several (n = 32) maternal-zygotic (MZ)dnah10 ut28/ut28 ; Tg(Foxj1a:Arl13B-GFP ut22/+ ) mutant larvae (generated by dnah10 ut28/ut28 x dnah10 ut28/ut28 ; Tg(Foxj1a:Arl13B-GFP ut22/+ mutant cross) at de ned regions of the ventricular system including: (i) the hindbrain ventricle, and within the central canal both at (ii) the vent/cloaca; and (iii) at the tip of the tail (Fig. 4K). When assessed at 3dpf, the majority of these mutant zebra sh exhibited dysfunctional cilia motility in the cerebral ventricle portion of the central canal (n = 23; 72%) ( Fig. 4K), while cilia motility within the central canal, at vent or in the tail, were less affected in aggregate at this time. To evaluate whether the dnah10 ut28/ut28 mutation affected overall cilia beat frequency, we sorted high-speed confocal imaging (33 Hz) of MZdnah10 ut28/ut28 ; Tg(Foxj1a:Arl13B-GFP) ut22/+ into two categories, "immobile cilia" and "motile cilia". As expected, the analysis of beat frequency of "immobile cilia" was zero (Fig. 4L, S4B). Similar analysis of beat frequency of the seemingly normal "motile cilia" in dnah10 ut28/ut28 mutants was reduced (6.2 ± 1.6 Hz, n = 22) compared with WT transgenic zebra sh (7.4 ± 2.3 Hz, n = 27), however this reduction was not signi cant at 3 dpf (Fig. 4L). Similar analysis of the beat frequency of motile cilia of WT Tg(Foxj1a:Arl13B-GFP) ut22/+ transgenic zebra sh found no signi cant differences in cilia beat frequency when comparing all three different locations of the hindbrain ventricle and within the central canal at vent and tail tip (n = 9; each region) ( Figure S4D).
After imaging, we reared all individuals (n = 32) for longitudinal analysis of the onset body axis curvatures during larval development. We observed that the majority of MZdnah10 ut28/ut28 ; Tg(Foxj1a:Arl13B-GFP) ut22/+ larvae displayed body curvature and scoliosis phenotypes by 28 dpf (90.4%; n = 19) (Fig. 4M), which was similar to our observations of incomplete penetrance of scoliosis observed in MZdnah10 ut28/ut28 mutants from independent clutches (Fig. 3G). Altogether, these results demonstrated that dnah10 is important for the motility of motile cilia in the brain and central canal, albeit expressed in stochastically in speci c regions at 3dpf and that this loss of central canal motile cilia beating is associated with the onset and progression of body curvature and scoliosis phenotypes in zebra sh.
dnah10 is required for maintaining Reissner Fiber assembly The Reissner ber (RF) has recently been shown to play an important role morphogenesis of a straight body axis during embryonic development (Cantaut-Belarif, Sternberg et al. 2018) and its continued assembly during larval development is essential for the maintenance of a straight spine in zebra sh (Rose, Pompili et al. 2020, Troutwine, Gontarz et al. 2020. To test whether motile ciliary defects in dnah10 ut2/ut28 mutant zebra sh are also associated with RF disassembly, we rst utilized immuno uorescence using antiserum raised against bovine RF (AFRU) in dnah10 ut28/ut28 mutant and WT zebra sh at 1, 5, and 10 dpf. As expected, WT zebra sh displayed the formation of a straight RF at all developmental timepoints (Figs. 5A-C). In contrast, the organization of the RF in dnah10 ut28/ut28 mutants at 1 dpf was already observed to be defective (Fig. 5A). The diffuse and discontinuous AFRU staining in dnah10 ut28/ut28 mutants at 1 dpf shows that the Reissner material is synthesized and secreted from the Sub-Commissural Organ (SCO), but is able to form a tightly polymerized RF within the central canal. By 5 and 10 dpf, we consistently observed more severe RF disassembly or the complete absence of a RF in dnah10 ut28/ut28 mutants, in contrast to the obvious RF structure in the WT animals ( Fig. 5B and C). Categorical analysis of the RF disassembly demonstrated a high degree of heterogeneity at 3 dpf (n = 25) including: (i) the formation of Reissner material as a bolus or punctate (n = 9, 36%); (ii) an intact RF which appeared distorted (n = 3, 12%); (iii) disassembly of the ber but with diffuse RF material lling the central canal (n = 8, 32%); or a complete absence of RF staining in the central canal (n = 5, 20%) ( Fig. 5D and E). These results, together with our analysis of cilia motility, demonstrate that the defects in RF correlated with the onset of cilia dysfunction in dnah10 ut28/ut28 mutants and suggests that the heterogeneity of RF disassembly may be related to the heterogeneity of cilia motility defects outlined above (Fig. 4K).
To facilitate the visualization of RF dynamics in dnah10 ut28/ut28 mutants we crossed these to the scospondin-GFP ut24 gene fusion knock-in zebra sh line (Troutwine, Gontarz et al. 2020). scospondin-GFP ut24 embryos display a uorescently labeled RF, which continually moves in a rostral-caudal direction (Fig. 6A, A'; Video S4). In contrast, the RF moved in a disorganized manner in scospondin-GFP ut24 ;dnah10 ut28/ut28 mutant embryos (Fig. 6B, C; Video S5). Analogous to RF defects demonstrated in xed, AFRU immuno-stained dnah10 ut28/ut28 mutant animals (Fig. 5), we observed heterogenous RF disassembly patterns scospondin-GFP ut24 ;dnah10 ut28/ut28 mutants in vivo (Video S6 and Fig. 6H) including: (i) the formation of bolus or punctate Reissner material ling the central canal (n = 5, 17%) (Fig. 6D); (ii) the complete absence of ber or Reissner material in the central canal (n = 4, 14%) (Fig. 6E); (iii) absence of a ber with only diffuse Reissner material lling the central canal (n = 14, 48%) (Fig. 6F); or (iv) a continuous RF surrounded with diffuse Reissner material (n = 6, 21%) (Fig. 6G). We observed heterogeneity in the ability of the RF to move along the central canal including moving in non-linear, meandering patterns (Fig. 6J), or slowly traversing in a rostral direction, or simply not moving within the central canal (Fig. 6L, M and Video S7). Altogether, these data demonstrate that alterations in the assembly of the RF and it typical rostral-caudal dynamics are caused by the loss of dnah10-dependent cilia motility lining the brain and central canal.

Discussion
In this study, we report exome sequencing results from 195 AIS patients. We identi ed a signi cant association of variants in axonemal dynein heavy chain genes in AIS patients compared to controls. An exome sequencing study in multigenerational families with idiopathic scoliosis also implicated components of the microtubule cytoskeleton including DNAH6 and DNAH8 (Baschal, Terhune et al. 2018). We generated a stable dnah10 ut28/ut28 mutant zebra sh line, which models aspects of AIS including late-onset, atypical three-dimensional spine curvatures without patterning defects. We go on to demonstrate that dnah10 ut28/ut28 mutants display heterogenous defects in cilia motility throughout the brain and spinal canal, which is consistent with heterogenous RF disassembly within the central canal and with the onset of body curvatures and adult-viable scoliosis. Altogether these results suggest that mutations in axonemal dynein genes may be involved in AIS susceptibility in humans.
Recent WES and GWAS have identi ed several candidate genes associated with AIS. However, only a few of these candidate genes have been veri ed in animal models. Of the 28 reported AIS candidate genes previously identi ed in genetics studies, only four genes (GPR126, MAPK7, POC5 and SLC39A8) have provided robust validation in an animal model (Karner, Long et al. 2015, Patten, Margaritte-Jeannin et al. 2015, Gao, Chen et al. 2017, Haller, McCall et al. 2018. Forward genetic screens in zebra sh have identi ed numerous scoliosis phenotype-causing genes, but few of these genes correspond to candidate genes implicated in human AIS Ciruna 2017, Liu and. Previous studies have also demonstrated that defects in ependymal cell motile cilia formation or beating after disruption of ciliary components (i.e. ptk7a, kif6, cfap298) (Buchan, Gray et al. 2014, Hayes, Gao et al. 2014, Grimes, Boswell et al. 2016) and the subsequent disassembly of the RF as in scospondin mutants (Rose, Pompili et al. 2020, Troutwine, Gontarz et al. 2020 are key initiating factors causing scoliosis in zebra sh. Here, we extend this model demonstrating an association of several axonemal dynein genes in human AIS and that disruption of one of these genes, dnah10, leads to scoliosis in zebra sh. Thus supporting a model where disruption epndeymal cell cilia beating in the hindbrain and spinal canal may underly the pathogenesis of human AIS. The most signi cantly enriched gene family of the 256 genes ltered using the gene-burden test was the axonemal dynein family. Axonemal dyneins are components of the inner and outer arms of motile cilia (Kamiya and Yagi 2014). In humans, variants in DNAH1, DNAH9, DNAH11, DNAH17 have been implicated in ciliary dyskinesia and spermatogenic failure (www.omim.org). DNAH10 has been associated with heterotaxy syndrome a congenital disorder characterized by malformation of left-right body asymmetry and associated with defective motile nodal cilia (Grimes and Burdine 2017). A previous study showed that disruption of dnah10 by an early terminating truncating mutation caused heart looping defects and randomization of left-right patterning genes (Liu, Cao et al. 2018). In contrast, we observed no obvious heterotaxy phenotypes in dnah10 ut28/ut28 mutant zebra sh reported here, most likely due to the hypomorphic nature of the C-terminal truncating mutation we generated, which could allow for residual protein function during early embryogenesis. In support of this model, we also consistently observed incomplete penetrance of scoliosis, and heterogenous defects in cilia motility and RF disassembly in MZdnah10 ut28/ut28 mutant animals where both the egg and sperm should only have mutant forms of the Dnah10 protein during embryogenesis. In addition, during the mulitplex CRISPR/Cas9 screening of the dnah10 locus using four guide RNAs we consistently observed animals displaying severe curled-tail down, larval lethal phenotypes, which is commonly associated with cilia defects in zebra sh. While, in the same injected clutch we also recovered milder body curvature and adult-viable scoliosis phenotypes associated with the stable dnah10 ut28/ut28 mutants reported here. A similar variance in the severity of body curvature phenotypes was observed between scospondin hypomorphic (Troutwine, Gontarz et al. 2020) and null (Cantaut-Belarif, Sternberg et al. 2018) mutant zebra sh and again in the comparison of the loss of maternal versus zygotic Ptk7a function (Hayes, Gao et al. 2014) are established in zebra sh. Altogther, this demonstrates that hypomorphic mutations in otherwise essential genes can display a range of body curvature phenotypes in zebra sh, which should inform the characterization of future human genetics studies of AIS.
Indeed, seminal human genetics studies have demonstrated that increased polygenic burden of rare variants across extracellular matrix genes in human AIS (Haller, Alvarado et al. 2016). Our current analysis of the Cterminal truncating dnah10 ut28/ut28 mutant zebra sh may suggest that milder mutations in axonemal dynein gene contribute to increase susceptibility of AIS, in contrast more deleterious mutations in the same gene may cause more severe motile cilia disorders attributed to primary cilia dyskinesia in humans and animal models.
Interstingly, the incidence of scoliosis is higher in primary cilia dyskinesia patients (Knowles, Daniels et al. 2013) and there is a strong association between organ symmetry and the direction of the spine curvatures reported in heterotaxic situs inversus patients, suggesting that left-right patterning has a key role in the regulation on spine morphogenesis (Schlösser, Semple et al. 2017). Further testing of this model of hypomorphic mutations underlying AIS by targeted genome editing of candidate human alleles of axonemal dynein genes associated with AIS is warranted.
We demonstrate that deleterious variants even with mild effect size in axonemal dynein genes may contribute to AIS susceptibility. However, replication of these variants in larger, multi-ethnic cohorts is necessary to help elucidate the inheritance pattern of these variants and to identify other genetic interactions and susceptibility loci for AIS. In addition, to axonemal dynein genes we also found by KEGG pathway analysis that ECM-receptor interaction related genes, which included three collagen genes (COL6A2, COL6A3, COL6A5, and COL6A6) were associated in our patient cohort. Another genetic study using WES performed in a European descent AIS cohort found that rare variants in COL6A3 were enriched in 435 AIS cases (Haller, Alvarado et al. 2016). We also identi ed two previously reported AIS candidate genes, HSPG2 and CELSR2, among our 256 candidate genes (Baschal, Wethey et al. 2014, Einarsdottir, Grauers et al. 2017).
In conclusion, our study has identi ed a signi cant accumulation of variants in axonemal dynein genes in Han Chinese AIS patients (Table S11) and showed that a C-terminal truncating mutation of dnah10 in zebra sh led to an adult-viable scoliosis phenotype by affecting motile cilia physiology and RF assembly. Further studies are needed to determine whether other components of motile cilia physiology may be involved in susceptibility of AIS in humans.