ABSTRACT
Among autistic individuals, a subphenotype with brain enlargement disproportionate to height (autism with disproportionate megalencephaly, or ASD-DM) seen at three years of age is associated with co-occurring intellectual disability and poorer prognoses later in life. However, many of the genes contributing to ASD-DM have yet to be delineated. In this study, we aim to identify additional ASD-DM associated genes to better define the genetic etiology of this subphenotype of autism. Here, we expand the previously studied sample size of ASD-DM individuals ten-fold by including probands from the Autism Phenome Project and Simons Simplex Collection, totaling 766 autistic individuals meeting the criteria for megalencephaly or macrocephaly and revealing 153 candidate ASD-DM genes harboring de novo protein-impacting variants. Our findings include thirteen high confidence autism genes and seven genes previously associated with DM. Five impacted genes have previously been associated with both autism and DM, including CHD8 and PTEN. By performing functional network analysis, we also narrowed in on additional candidate genes, including one previously implicated in ASD-DM (PIK3CA) as well as 184 additional genes previously implicated in ASD or DM alone. Using zebrafish as a model, we performed CRISPR gene editing to generate knockout animals for seven of candidate genes and assessed head-size and induced seizure activity differences. From this analysis, we identified significant morphological changes in zebrafish knockout models of two genes, ythdf2 and ryr3. While zebrafish knockout mutants model haploinsufficiency of assayed genes, we identified a de novo tandem duplication impacting YTHDF2 in an ASD-DM proband. Therefore, we also transiently overexpressed YTHDF2 by injection of in vitro transcribed human mRNA to simulate the patient-identified duplication. Following this, we observed increased head and brain size in the YTHDF2 overexpression zebrafish matching that of the proband, providing a functional link between YTHDF2 mutations and DM. Though discovery of additional mutations of ASD-DM candidate genes are required in order to fully elucidate the genetics associated with this severe form of autism, our study was able to provide support for YTHDF2 as a novel putative ASD-DM gene.
INTRODUCTION
Autism is a group of neurodevelopmental traits characterized by difficulties with communication, social interactions, and behavioral challenges, prevalent in 1 out of 54 children in the United States 1,2. Autism is highly heritable, with 50–90% of cases estimated to be driven by genetics alone 3–5. Autism is also highly heterogeneous with recent large-scale whole exome sequencing (WES) of >63,000 autistic probands identifying 125 high confidence autism genes, with the predicted number of genes left to be discovered exceeding 1,000 6–8. In particular, leveraging genomic data from autism families—including parents and unaffected siblings—in the Simons Simplex Collection (SSC) 9 has identified coding de novo variants estimated to contribute to 30% of diagnoses 10–12. More recently, whole genome sequencing (WGS) of SSC de novo noncoding mutations implicates risk in an additional 4.3% of autism cases 13,14. Despite the combined efforts to sequence tens of thousands of genomes, known genes still only account for 10-20% of cases, and further work is required to fully elucidate genes and pathways contributing to autism etiology 15. Combining de novo variation with autism sub-phenotyping has been used to address the heterogeneity of autism and identify susceptibility loci for comorbid phenotypes in an acute way 16.
Brain enlargement that is disproportionate to height, known as disproportionate megalencephaly (DM), is enriched in autistic probands with 15% of autistic boys falling under the DM subphenotype (ASD-DM) compared to 6% in typically developing (TD) boys 17. This comorbidity is associated with more severe cognitive phenotypes, including lower IQ and language use, as well as higher rates of language regression 18–20. This robust enrichment and distinct presentation support DM as a sub-phenotype of ASD, likely due to a shared genetic etiology between autism and DM. While a handful of genes have been associated with DM—including known autism genes impacting cell cycle and proliferation during embryonic development (e.g., CHD8 and PTEN)—mutations of known candidate genes make up only 3% of megalencephaly in autism probands. This leaves the genetic etiology of a majority of ASD-DM cases undiscovered 21–24. A recent study using WES from 46 autistic families with macrocephaly (ASD-M)—defined as >2 standard deviations above the mean in head circumference for TD sex and age-matched children—successfully identified mutations in one novel and several known autism candidate genes 25, demonstrating the power of sub-phenotyping ASD-DM leading to genetic discoveries even for reduced sample size.
Zebrafish (Danio rerio) are an attractive model for studying neurodevelopmental traits given their rapid development, large number of progeny, transparent bodies, and that ∼70% of gene orthologs are shared with humans 26,27. Previous studies of known ASD-DM genes recapitulate macrocephaly and DM phenotypes in zebrafish knockdown and knockout experiments for CHD8 and KMT2E, respectively 28,29. This method of knocking down candidate ASD-M genes has also been used systematically to identify the contributing gene in the chromosome 16p11.2 locus in zebrafish 30. Further, novel technologies such as the VAST BioImaging System allow for the rapid characterization of zebrafish knockout models through the generation of high resolution standardized images 31. We recently demonstrated the utility of this approach by examining the knockdown of two genes associated with autism and microcephaly, SLC7A5 and SYNGAP1, by assessing head-size phenotypes in CRISPR-generated zebrafish knockout line embryos at three and five days post fertilization (dpf) 32.
In this study, we leveraged high-coverage WGS data from two cohorts, 11 ASD-DM probands from the UC Davis MIND Institute Autism Phenome Project (APP) specifically identified using magnetic resonance imaging (MRI) data at around three years of age, and 755 ASD-M probands with head circumference data available from the SSC cohort. Together, this represents a >10-fold increase in probands compared to the previous largest study of increased head circumference associated genes in ASD 25. Using this sub-phenotype-to-genotype analytic strategy, we identified candidate ASD-DM and ASD-M genes harboring de novo loss-of-function (LoF) variants, including CHD8 and PTEN, and subsequently functionally validated a subset of candidate genes using zebrafish. From this, we narrowed in on two genes impacting larval head size. One of these genes, YTHDF2, also lead to macrocephaly in zebrafish following embryonic microinjection of mRNA, recapitulating the phenotype seen in the ASD-DM proband harboring a partial tandem duplication of YTHDF2. Together our sub-phenotyping approach provides a powerful strategy to identify novel ASD-DM candidate genes and validate their role in brain development in a zebrafish model system.
RESULTS
ASD-DM Candidate Gene Discovery
ASD-DM individuals were recruited via the UC Davis APP—a longitudinal study focused on the identification of ASD-subphenotypes 17,33,34. Using MRI data from the study entry time point (2–31/2 years of age) 33, we identified 11/89 WGS sequenced individuals in the APP cohort that met the criteria for ASD-DM, defined as a cerebral volume to height ratio >1.5 standard deviations above the mean compared to TD age-matched controls. Through a collaboration with MSSNG 35,36, WGS and variant identification/annotation was performed for the autistic probands and a subset of family members, for which we also had blood specimens, including six trios and five non-trio probands yielding over 200 thousand variants. From this, we identified two exonic, de novo, LoF variants from trio families, including one single-nucleotide variant (SNV) splice-site variant impacting RYR3 and one 109-kbp duplication of YTHDF2 and GMEB1. An additional ten variants were found in non-trio proband data to be exonic, LoF, and rare (not previously recorded in dbSNP) 37.
We also examined variants from the ASD-DM probands where we do not have sequenced parental genomes to determine de novo variant status, but may harbor rare variants associated with their ASD-DM phenotype. From these five individuals, we identified a proband harboring a chromosome 1q21.1 microduplication, a CNV previously associated with ASD-DM 38, and a single proband with variants in CHD8 and KMT2E 39. We also uncovered ten ultra-rare variants in these probands, defined as previously undescribed in dbSNP (dbSNP). Of these, three impacted genes have SFARI scores of 3S or above (KMT2E, RPS6KA5, and TTN), an additional three genes have known neuronal functions (DMBT1, IARS2, FGF12), and one gene was found recurrently carrying variants in two probands (SPANXN4) 40.
The SSC resource consists of genetic and phenotypic information from trios and quads of simplex autism families. SSC head circumference and age data was used to determine ASD-M status (head circumference >1.5 standard deviations above the mean for typically developing sex and age-matched children) for 755 of 1,846 SSC probands (40%) 9 (Supplementary Table 1) (Supplementary Figure 1). De novo variants from WGS of ASD-M trios and quads were assembled from previously published SSC de novo exonic variant calls 13. In total, 157 ASD-M de novo variants impacted coding exons and were predicted to be LoF, overlapping a total of 150 genes (Figure 1) (Supplementary Table 2). Of note, five genes were found to contain variants in multiple probands, including GALNT18, KDM6B, LTN1, RERE, and WDFY3, as well as CHD8, which was disrupted in three distinct probands. For SSC probands, only SNVs and indels were considered. Due to the lack of MRI data for SSC participants, for this study we used ASD-M as a proxy for ASD-DM.
A. A histogram representing the number of SSC probands v. head circumference percentiles shows a skew towards larger head-sizes compared to age and sex-matched typically developing children. The red bar designates those meeting the criteria for macrocephaly. The dashed line represents the distribution mean. B. ASD-DM and ASD-M genes listed by their identified proband’s head circumference percentiles show genes previously associated with ASD-DM (first grey quadrant) are more likely to be associated with a higher head circumference percentiles than genes previously associated with autism (second white quadrant) and DM (third grey quadrant) alone. Color represents the macrocephaly type. DM, disproportionate macrocephaly; RM, relative macrocephaly; SO, somatic overgrowth.
In total, we identified 153 genes containing a LoF variant in the APP ASD-DM and SSC ASD-M datasets. Rates of harboring a LoF de novo variant in ASD-DM probands were in line with previous predictions (19.4%) and nominally enriched compared to TD SSC siblings (16.5%), though not statistically significant (chi-squared p-value = 0.1) 10,25. Over a third of identified candidate genes (53/153) had a pLI score of > 0.9, suggesting they are intolerant to variation 41 (Supplementary Table 2). Examining de novo missense variants, which previously exhibit overall enrichments in affected probands versus unaffected siblings 42,43, we did not observe an enrichment in our LoF candidate genes in ASD-M probands compared to TD siblings of ASD-M probands, ASD-N probands, and siblings of ASD-N probands (Student’s T-test, p-values = 0.16, 0.29, 0.76).
ASD-DM Candidate Gene Discovery Using Network Analysis
Identifying shared patterns of molecular functions and ontologies of impacted ASD-DM genes may point to additional gene candidates. Due to the highly heterogeneous nature of ASD, this type of analysis expands our ability to identify disrupted biological mechanisms and spatio-temporal expression patterns implicated in autism 24. Here, we used as seeds 166 previously-known and identified-in-this-study ASD-DM genes to identify active interactions using the STRING database (Figure 2) 25,44. This analysis uncovered ontology groups enriched in our dataset previously reported for ASD, including proteins involved in histone modification and chromatin organization, transcription factors, cell signaling (e.g., SMAD and E-box binding), functions key to neuronal activity (e.g., sodium and calcium ion transport and glutamate receptor binding), cell adhesion and cytoskeletal proteins, and mRNA binding (Supplementary Table 3) 45–48. Out of our original ASD-DM candidate seed genes 41.5% (69/166) fall under one of these ontologies.
ASD-DM candidate genes from SSC (teal), APP (purple), and Wu (navy) probands are connected in a network via active interactions as determined by STRING 44. Background colors represent shared gene ontology (GO) molecular functions. Disconnected gene nodes are not included.
In order to identify ontologies specific to ASD-DM compared to ASD and DM alone we used the Database for Annotations, Visualization, and Integrated Discovery (DAVID) to identify unique ontologies enriched by the 166 known and candidate ASD-DM genes compared to ontologies enriched in SSC ASD-N proband LoF genes, and genes previously associated with DM 49,50. ASD-DM is uniquely enriched for terms such as intellectual disability, synapse assembly and long term synaptic depression, histone methyltransferase activity, cytoskeletal structure (spectrin repeats) (Supplementary Table 3). Though not reaching statistical significance, ASD-DM is also weakly enriched for mRNA binding with 2/19 previously annotated m6A-mRNA regulating genes found in ASD-DM (YTHDC1, YTHDF2) with no hits in ASD-N or in TD SSC siblings 51. ASD-DM was also not enriched for genes in the MTOR pathway (NF1 and PTEN) compared to ASD-N (GRB10, PPP2R5D, RICTOR, and TSC1) 52.
We next sought to identify putative additional candidate ASD-DM by expanding our network to include the top ten interactors for each ASD-DM seed gene (Supplementary Table 4). This list of top ten interactors included genes defined by STRING as having known protein interactions, shared homology, and co-expression patterns 44. Of the top ten ASD-DM candidate gene interactors, 28% (518/1826) also fall under one of the ontologies found in our ASD-DM network. Interestingly, one of these genes has previously been associated with both autism and DM individually, PIK3CA (Supplementary Table 5). PIK3CA is a catalytic subunit of the mTOR pathway previously found to be associated with developmental delay and DM, including one diagnosed with autism 53.
In this top ten interactor set, twenty genes are high-confidence autism genes not previously associated with DM (ANK2, ASXL3, CTNNB1, CUL3, DLG4, DYRK1A, GNAI1, GRIN2B, KCNMA1, KMT2A, NCOA1, NIPBL, NRXN1, PHF12, POGZ, PPP1R9B, SIN3A, SMARCC2, TBL1XR1, UBR1). Eighteen additional genes from this interactor set have been implicated in DM and implicated as SFARI putative autism candidate genes (ANK3, CHD2, CHD3, FRMPD4, HCFC1, HDAC4, HRAS, HUWE1, PAK1, PIK3R2, RAC1, SETD1A, SLC25A1, SMAD4, TBL1X, TRIO, USP7, and USP9X), and 184 more have been implicated in DM or have a SFARI score. Especially promising are the 21 genes that contain missense variants in the SSC ASD-M probands, but do not contain de novo missense variants in ASD-N probands or their TD siblings, including ABI2, ANK3, SRC, SRCAP, ATP12A, BAIAP2, CHD13, CH815, FGG, JUP, KDM2A, KIF20A, MAPK8, PDGFRB, RING1, SCN4A, SHANK1, SMC3, TCF3, WDR5, ZC3H3. Together, network analysis and ontology point to these genes as promising ASD-DM candidate genes going forward.
Validation of Candidate Genes in a Zebrafish Model System
In order to determine if our candidate ASD-DM genes were associated with a head-size phenotype, we tested seven by CRISPR knockout in zebrafish (Supplementary Figures 2 & 3). We focused our analysis on the previously unstudied candidates from the APP cohort, including the three genes impacted by de novo variants, RYR3, GMEB1, and YTHDF2. It is known that ∼20% of zebrafish genes conserved with human are also duplicated; we therefore prioritized candidate ASD-DM genes harboring rare variants in our ASD-DM non-trio dataset that have a single ortholog in zebrafish (IARS2 and RPS6KA) 54. Additionally, we included two genes from the SSC ASD-M cohort, CHD8, in which CRISPR F0 mosaic knockout has previously been shown to lead to increased interorbital distance in zebrafish embryos 28, and FAM91A1, a gene previously not implicated in ASD. Finally, as controls, we included three genes found in SSC ASD-N probands for which we do not expect to see a head-size phenotype—HEPACAM2, PAX5, and SCP2A.
We generated CRISPR knockout F0 embryos by microinjection of four guide RNAs (gRNAs) targeting exonic regions for each assessed gene previously shown to result in near complete mosaic knockout of target genes with little off-target effects 25,32,55. We then used the VAST BioImaging System paired with FishInspector software to perform automated measurements of embryo area, length, and distance between the center of the eyes (as a proxy for head size), at 3 days post fertilization (dpf). From this analysis, we identified two genes (ryr3 and ythdf2) showing significant head-size differences in zebrafish knockout embryos compared to negative scrambled injection controls (Wilcoxon T-test p-values < 0.01) (Figure 3). In both cases, knockout embryos exhibited reduced head size consistent with microcephaly, with ryr3 showing a mean reduction of 3.3% and ythdf2 showing a mean reduction of 3.8%. RYR3 encodes for a ryanodine receptor responsible for calcium transport in muscle and brain tissue, and YTHDF2 encodes for part of the m6A-mRNA degradation complex 56,57. Though ythdf2 also displayed a reduction of body area and length, this size difference is not due to overall developmental delay measured by head trunk angle (Supplementary Figure 4). Interestingly, only genes harboring de novo variants from the ASD-DM APP cohort showed head-size phenotypes, though we did not observe a head-size phenotype in our chd8 knockout embryos using our morphometric measurements, counter to previously published results 28.
Phenotyping knockout zebrafish morphometric measurements using the VAST BioImaging System, we identified ryr3 and ythdf2 as having reduced head size compared to a negative scrambled control. Assessing for increased seizure-like activity in the presence of PTZ using ZebraBox technology we identified three genes associated with increased high speed events, and 3 associated with decreased high speed events. Circle size represents the inverse Wilcoxon T-test p-value. Circle color represents the fold change compared to a negative scrambled control.
As increased prevalence of seizures is highly enriched in autism 58, we assessed our knockout embryos for increased susceptibility to drug-induced seizures 59 by treating them with GABA antagonist pentylenetetrazole (PTZ; 5 mM) at 5 dpf 32,60. High speed events, corresponding to seizure-like movements were recorded using ZebraBox 61 and quantified using published software 62 (Figure 3). This analysis revealed three gene knockout models associated with an increase in high speed events versus scrambled controls, fam91a1, gmeb1, and hepacam2. Of note, the human chromosome 7q21.3 microdeletion syndrome containing HEPACAM2, an immunoglobulin gene responsible for cell adhesion, which showed the highest fold increase in high speed events compared to controls in our zebrafish assay, have previously been associated with myoclonic seizures in humans 25,63. Mutations in FAM91A1, responsible for Golgi protein trafficking, have been linked to dysregulated electrophysiological brain activity 64,65. Though previously unassociated with seizures, GMEB1 - a caspase activation and apoptosis inhibitor, does have a known association to schizophrenia, which shares a historical link with epilepsy 66–68. Moving forward, we focused specifically on YTHDF2, found to be nearly fully duplicated in an APP ASD-DM proband and leading to reduced head size in knockout zebrafish embryos. This due to its clear mechanism of action - gain-of-function due to overexpression.
YTHDF2 in ASD-DM
Based on the findings that CRISPR loss-of-function leads to microcephaly for ythdf2 and increased seizure activity for gmeb1, we sought to better characterize the de novo 109 kbp duplication identified in an APP proband with ASD-DM. We validated this CNV using sequence read depth (QuicK-mer2) 69,70 (Figure 4A). Split reads falling at the identified breakpoints of this copy-number variant (CNV) suggested a tandem duplication harboring the entire GMEB1 gene and the first five of six exons of YTHDF2 inserted at the 3’ untranslated region (UTR) of the noncoding divergent transcript of TAF12, directly upstream of GMEB1 (Figure 4B and C).
A. QuickMer2 plot showing the predicted copy number of the APP proband harboring a de novo duplication a 110kb duplication on chromosome 1 compared to their parents harboring only 2 copies. B. IGV plot showing discordant reads in the APP proband supporting a tandem duplication. C. An illustration of the 110kb tandem duplication on chromosome 1 in an APP proband encompassing GMEB1 and all but the last exon of YTHDF2. D. Head size of mRNA injected zebrafish, estimated via telencephalon measurement, are significantly increased after YTHDF2 mRNA microinjection into zebrafish embryos (Wilcoxon T-test p-value = 0.0013), and unchanged in GMEB1 mRNA injected embryos at 3dpf (Wilcoxon T-test p-value = 0.8). E. Knockout and mRNA injected zebrafish harboring a pan neuronal marker reveal brain size differences at 3dpf. Knockout embryos show significantly decreased midbrain volume (Wilcoxon T-test, p-value = 0.003). mRNA injected embryos show both significantly increased midbrain (Wilcoxon T-test, p-value = 0.007) and forebrain (Wilcoxon T-test, p-value = 0.0005). F. Representative control, knockout, and mRNA injected zebrafish images from wildtype embryos imaged with the VAST bioimaging system and HuC embryos harboring a pan neuronal fluorescent marker. Scale bar is 100 uM.
Using available microarray data produced from mRNA derived from whole venous blood, we found that both GMEB1 and YTHDF2 exhibited significantly increased expression in the APP proband harboring the duplication compared to other APP participants (Student’s T-test p-value < 0.0001) 71. GMEB1 is an auxiliary factor in parvovirus replication known to inhibit apoptosis in neurons and previously associated with schizophrenia 67,72 and YTHDF2 is a member of the m6A-containing mRNA degradation complex known to be downregulated in neuronal fate determination 57, making both of these attractive potential ASD-DM candidate genes.
We hypothesized that the knockout of ythdf2 leading to microcephaly in zebrafish larvae made its ortholog a top candidate for YTHDF2 overexpression leading to ASD-DM. Based on their known functions, duplication of either gene could plausibly result in neurodevelopmental effects. Therefore, we modeled increased expression of YTHDF2 and GMEB1 by microinjecting human in vitro transcribed mRNA into single-cell stage zebrafish embryos. We then assessed these ‘overexpression’ zebrafish embryos at 3 dpf using the VAST BioImaging System (Figure 4D). For YTHDF2, we observed increased telencephalon length compared to a dye injected negative control consistent with increased head size but no significant differences to body area or length, matching the proband phenotype. We did not observe any difference in head size for the GMEB1-injected larvae.
To verify the increased head size was a result of megalencephaly, we repeated the experiment in the zebrafish transgenic line HuC-eGFP, which harbor a GFP fluorescent pan-neuronal marker (Figure 4E and F). These embryos displayed brain size differences, with knockout embryos showing significantly reduced midbrain, and mRNA injected ‘overexpression’ embryos exhibiting significantly increased midbrain and forebrain compared to injection controls. Together, the knockout and mRNA injected zebrafish provide evidence that gain of YTHDF2 is associated with DM while its loss leads to microcephaly.
To better characterize ythdf2 neurodevelopmental functions, we sought to understand expression of this gene. Previous analyses of ythdf2 expression in zebrafish via in situ hybridization revealed widespread robust expression 73. This agrees with reported widespread expression in both humans and zebrafish, beginning early in development at 3 hpf in zebrafish embryos (Supplementary Figure 5) 74,75.
Not only is YTHDF2 highly conserved, with 95% gene conservation to mouse and 72% in zebrafish, it also has very few known annotated variants in humans, suggesting it plays a key role in development 76. YTHDF2 is predicted to be intolerant to mutation, with a pLI score of 1, and a LoF splice site variant found in the Gnomad Database found in a single male of Latino descent (GRCh38, chr1:28738331,C>G) 41 (Figure 5A). There are also no annotated CNVs in the 1000 Genomes Project or SSC overlapping YTHDF2, and no CNVs in ClinVar encompassing fewer than 16 genes (hg38, chr1:28493687-29242679; del), with some larger CNVs associated with probands diagnosed with intellectual disability and developmental delay 9,77–79. Notably, two individuals in the SSC database were found to harbor YTHDF SNVs, one female proband with a de novo missense variant (hg38, chr1:28742997,G>A), and one unaffected parent harboring a LoF splice site mutation (hg38, chr1:28737684,T>C), consistent with incomplete penetrance.
A. Known YTHDF2 variants are rare, and primarily found in autistic probands and family members. B. Proposed model of YTHDF2 deletion and duplication phenotypes. We hypothesize YTHDF2 deletion would lead to microcephaly due to lack of m6-mRNA degradation, extended cell cycle progression, and reduced neurogenesis. Inversely, YTHDF2 duplication would lead to megalencephaly following increased m6A-mRNA degradation, neural progenitor cell (NPC) overabundance, and increased neurogenesis.
DISCUSSION
The autism sub-phenotype of disproportionate megalencephaly (ASD-DM), which occurs in approximately 15% of autistic boys, is associated with lower language ability at age three and slower gains in IQ across early childhood resulting in a higher proportion with IQs in the range of intellectual disability by age six 17. Others have reported similar rates of macrocephaly in autism (ASD-M) 80. In this paper we have examined the genomes of 766 ASD-DM and ASD-M trios and quads from the APP and SSC cohorts to identify 153 ASD-DM candidate genes containing de novo LoF variants. We functionally tested seven of these candidate genes using a zebrafish CRISPR knockout model and discovered two of these genes to be associated with reduced head size in zebrafish: ryr3 and ythdf2. In the case of the ryr3 knockout, our findings of microcephaly are counter to those observed in the autistic proband with an identified LoF splice site variant of RYR3 resulting in DM (chr15:33634585; G>T). Further experimentation delineating the molecular mechanisms of this proband variant will be necessary in order to determine if reduced head size in knockout zebrafish embryos is due to a potential gain-of-function mechanism or intra-species differences of ryr3 function. For ythdf2, to appropriately model the identified patient YTHDF2 duplication for which we hypothesize gene gain-of-function effects, we overexpressed YTHDF2 in zebrafish, recapitulating the increased head-size phenotype. Together, this paper introduces two novel ASD-DM candidate genes associated with head-size phenotypes in a zebrafish model system, 141 novel unvalidated ASD-DM candidate genes, and establishes a rapid pipeline for the identification and validation of these genes. In addition to YTHDF2 and RYR3, the candidate genes identified in this paper as harboring de novo LoF variation in ASD-DM APP and ASD-M SSC probands also make attractive targets for future phenotypic validation. This list can be further prioritized by choosing to follow-up on genes with similar functions to known ontologies in autism and ASD-DM, as well as ranking head-size phenotypes in the probands. With this expanded list of ASD-DM candidate genes, network analysis can now be leveraged to identify additional candidate genes with similar gene functions to known ASD-DM genes.
We hypothesize that the strong enrichment between autism and DM point to a shared genetic etiology between these two traits. As 43% of de novo likely-gene-disrupting variants in probands have been estimated to contribute to an autism diagnosis, we expect many genes identified in this study to contribute to ASD-DM or autism alone 10. Out of the 153 genes we identified harboring a variant in an ASD-DM or ASD-M proband, five have previously been associated with both autism and DM (Table 1). These include well-characterized ASD-DM genes CHD8 and PTEN. CHD8, a chromatin remodeler important in early brain development for which we recovered three distinct variants, is one of the genes most strongly associated with autism and intellectual disability and 28,81. PTEN, another high confidence autism gene, is a tumor suppressor gene known to function in cell proliferation 82. To date, >26 variants of CHD8 and >40 variants of PTEN have been identified in autistic probands, and both genes have been estimated to have an 85% penetrance of overgrowth and macrocephaly phenotypes 83,84. Together, CHD8 and PTEN variants have been estimated to contribute to up to 15% of all ASD-M cases 25. Additionally, we recovered rarer high confidence autism genes also associated with macrocephaly, KMT2E, SHANK3, and WDFY3 6,50. Over 31 variants have been found in KMT2E, which plays a role in regulating histone methylation, in cases of ASD, with an estimated 45% macrocephaly penetrance 85. Mutations in SHANK3, a postsynaptic scaffolding protein, are known to be causal in the ASD-M syndromic condition Phelan–McDermid syndrome 86. WDFY3, which we recovered twice in our cohorts, encodes for an autophagy scaffolding protein 87. Interestingly, LoF variants in WDFY3 are associated with macrocephaly, while a likely gain-of-function WDFY3 variant has been found in an autistic individual with microcephaly.
We recovered thirteen additional high confidence autism genes not previously associated with DM 6,50 (Table 1). For 5/13 of these high confidence autism genes, including CELF4, CTCF, KDM5B, KDM6B, and TRIP12, we identified no LoF variants in the ASD-N probands and a variant has been previously identified in at least one ASD-DM proband, suggesting a shared etiology contributing to both autism and DM may exist for these genes but has yet to be robustly described 88–92. Of note, ARID1B was identified as a candidate ASD-M gene by a previous study 25. We identified an additional seven genes previously associated with DM, but not known to be a high confidence autism gene. Of these genes, 6/7, including LRP2, NFIB, PSMD12, RERE, SETD2, and WDFY3 have a SFARI score of two or above, suggesting they have been found in a case of ASD, but do not meet the criteria for a high confidence autism gene 40. Together, these 20 genes make attractive candidate genes, with shared gene functions that mirror high-confidence ASD-DM genes, such as cell signaling and neuronal processes (8/19) and histone and chromatin processes (5/20). Out of the remaining 128 genes not previously associated with autism or DM, 26 have SFARI scores two and above 40. Of these, we also identified 6 genes found to contain variants in more than 1 proband (Table 1). This includes the known ASD-DM chromatin remodeler CHD8 - found in 3 distinct probands, the autism associated demethylase KDM6B, the DM-associated retinoic acid signaling RERE and the autophagy-linked WDFY3, as well as two gene previously associated with either autism or DM include GALNT18 which functions in O-linked glycosylation and LTN1 which functions as a ubiquitin ligase 28,87,93–96.
In addition to the work described in this paper, previously published work supports roles of YTHDF2 in neurodevelopment. YTHDF2 is a m6A-specific RNA binding protein which designates m6A labeled RNA for degradation 97. It is known to bind over 3,000 transcripts, including targets previously associated with ASD, such as CREBBP 57. Loss of m6A-RNA degradation activity has been studied through the depletion of methyltransferase complex proteins, which leads to delayed neurogenesis and extended cell cycle progression of cortical neural progenitors—a mechanism associated with microcephaly 98,99. Additionally, the downregulation of YTHDF2 is required for neuronal fate determination and loss of YTHDF2 is associated with delayed mitotic entry 100,101. YTHDF2 is also known to function inversely to the ASD-syndrome gene FMRP, loss of which is associated with Fragile X Syndrome, which protects its target m6A-RNA from degradation by YTHDF2 102. Knockout mouse models of Ythdf2 have also shown decreased cortical thickness caused by decreased neurogenesis in early development 103. Meanwhile, studies in pigs have associated reduced Ythdf2 activity with inhibited porcine induced pluripotent stem cell pluripotency104. Together, this evidence suggests a model of YTHDF2 duplication in which neuronal fate determination is delayed, leading to the overabundance of neuronal progenitor stem cells followed by increased neurogenesis (Figure 5B).
Interestingly, we identified another gene coding for a m6A-specific RNA binding protein identified in the SSC ASD-M cohort, YTHDC1, an m6A-RNA reader that facilitates recruitment of splicing factors 105. Where gain of YTHDF2 may lead to a decrease in overall m6A-mRNA through increased degradation, conversely the loss of YTHDC1 would lead to a decrease of appropriate spliced m6A-mRNA 106. Further, viral connections between m6A-RNA regulators and brain size support the role of YTH domain protein family members in head size phenotypes. The m6A-labeled mRNA flavivirus ZIKA is associated with severe congenital microcephaly, with YTHDF2 found to bind and destabilize viral RNA 107,108. These antiviral functions are controlled in part by the METTL3 methyltransferase, which labels viral RNA for degradation, and whose knockout models are also associated with a reduced brain size in mice 109. Our network analysis also highlighted an enrichment of mRNA binding genes in general, including genes with known splicing activity such as HNRNPL, a mechanism previously found to be disrupted in autism 110. Together the m6A-mRNA pathway including YTH domain proteins and m6A methyltransferases and demethylases as well as mRNA binding genes will make for an interesting future area of study in regard to brain size phenotypes.
From the literature, twenty genes have been implicated as high confidence autism genes and have separately been implicated in DM 6,7,50. Together, the candidates identified in our analysis as well as previously identified ASD-M candidate genes cover 9/20 of these genes, suggesting further work will be required to capture the full genetic heterogeneity of ASD-DM 25.
In addition to an enrichment of both megalencephaly and macrocephaly in ASD, previous studies have found both to be correlated with decreased IQ and increased levels of language regression 17,111,112. Additionally, macrocephaly has previously been found to be highly correlated with megalencephaly, with an increased correlation in young children 113. Unfortunately, due to the dearth of MRI evidence, not all ASD-M probands included in this study will meet the criteria for ASD-DM. Though we expect the overlap to be significant, supported by the majority of APP ASD-DM probands also meeting the criteria for macrocephaly (82%), many of the candidate ASD-DM genes identified in this study may have arisen from false positive ASD-DM probands. Additionally, these genes may harbor de novo LoF by chance, with no contribution to a head or brain size phenotype in the proband. Going forward, it is paramount to validate these candidate ASD-DM genes as related to head or brain-size phenotypes. Additionally, we expect a fraction of cases to represent false negatives due to lack of sensitivity of the instrument, molecular processes such as genetic compensation, and inefficient gene knockout which may explain our lack of phenotype seen in CRISPR knockout chd8 fish compared to previous models 28,114. This result highlights a disparity in the CRISPR-VAST pipeline, suggesting other methods will be necessary to uncover false-negative ASD-DM genes from this analysis.
This research represents a significant increase in the number ASD-DM and ASD-M proband genomes analyzed in search of candidate genes. The 153 candidate ASD-DM genes introduced here greatly expand our knowledge of the genetic factors specifically contributing to this severe subphenotype of ASD. Additionally, our research establishes a roadmap for the rapid functional characterization of these putative ASD-DM, proving the methodology needed to validate these candidate genes going forward.
METHODS
Megalencephaly and Macrocephaly
APP probands and determinations of megalencephaly were previously determined as part of the APP study17. Acquisition of MRI data for megalencephaly measurements were made during natural nocturnal sleep for children at study enrollment (time point 1), between the ages of 2 and 3.5 115. Blood collected during this time point was sequenced using the to 30X coverage whole genome sequencing (WGS) through a collaboration with MSSNG 35,116. Raw data including FASTQ and VCF files can be accessed through the MSSNG access agreement: https://research.mss.ng.
Macrocephaly cases from the SSC were defined using a permissive cutoff of a head circumference >1.5 standard deviations (90%) above the mean of age matched controls 82. Age matched TD head circumference data 117 and height data 118 from ages 4-17 were derived from publicly available standards. For males in the SSC cohort 551/1601 (34%) met the criteria for macrocephaly, for females 108/245 (45%) met this criteria.
We identified three types of macrocephaly: (1) somatic overgrowth (SO) with head circumference and height percentiles >90%; (2) disproportionate macrocephaly (DM) with height percentiles over head circumference percentile < 0.7; and (3) relative macrocephaly (RM) with height percentiles over head circumference percentile > 0.7.
Variant Annotation
De novo variants were identified in APP probands as those unique to the proband and absent from either parent via string matching (grep -Fvxf). Variant mapping and type predictions were carried out as part of the MSSNG consortium pipeline 35,116. We considered LoF variants as those predicted to lead a frameshift, nonsense, or splice site mutation. De novo LoF variants were validated by IGV 119. Rare variants identified using dbSNP as those with a minor allele frequency (MAF) < 0.2% in all 5 1000 genome super populations 77,120. The presence of all rare and de novo variants identified in the APP cohort were validated by visual inspection of sequencing data via the Integrated Genomics Viewer (IGV) 70.
Network and Ontology Analyses
Network analysis of known ASD-DM genes and candidate ASD-DM genes from this study were completed using the STRING database and visualized via Cytoscape 44,121. Gene ontology analysis was completed for known ASD-DM genes and candidate ASD-DM genes as seed genes along with their top 10 gene interactors determined using the STRING database, similar to previous studies 24. Similar STRING database Molecular Function GO terminology were pooled. Gene ontology was completed using Database for Annotations, Visualization, and Integrated Discovery (DAVID) software 49. Human genes were used as background for ontology analyses.
Zebrafish CRISPR and mRNA models
gRNA were chosen using CRISPRScan for gRNA with a score of 35 or higher (Supplementary Table 6)122. In order to achieve 90% knockdown efficiency, we designed four gRNA per target gene as in Wu et al. 123. 2nmol crRNA were ordered from Integrated DNA Technologies (IDT). gRNA efficiency was tested via post-injection using a pooled genomic extraction of 4 embryos followed by 7.5% polyacrylamide gel visualization. Quantification of gRNA efficiency was carried out via amplicon sequencing and the CrispRvariants R package (Supplementary Figure 3, Supplementary Table 6) 124.
Injection mixes of ribonucleic protein (RNP) consisting of 4 pooled gRNA and Cas9 Nuclease, S pyogenes (New England Biolabs, M0386M) were prepared as previously described 32. Pooled gRNA were microinjected into single cell-stage NHGRI-1 zebrafish embryos to a volume of 0.5 nL/cell as previously described using a Pneumatic MPPI-2 Pressure Injector 125. Injection needles were made using Sutter Instrument Model P-97 Micropipette puller, with the settings Heat 680, Pull 90, Vel 90, Time 225, and thin wall glass capillaries (4m 1mm x 0.75 Fisher Scientific, 50-821-984). Scrambled injection RNP mix contained a single gRNA designed to have no target in the zebrafish genome.
Human mRNA was generated using YTHDF2 and GMEB1 cDNA plasmids (Horizon YTHDF2, MHS6278-202827242; GMEB1, MHS6278-202827172) 126 and prepared via the in vitro transcription kit mMESSAGE mMACHINE™ SP6 Transcription Kit (Thermo Fisher Scientific, AM1340). Injection mix of 100ng/uL mRNA and 0.05% phenol red were prepared as previously described and injected to a volume of 0.5 nL/cell127.
Zebrafish Morphometric Measurements
We took dorsal and ventral images of 3dpf embryos using the Union Biometrica VAST Bioimaging System with LP Sampler using the built-in camera and manufacturer settings 31. Zebrafish features were then identified and quantified from VAST using FishInspector software 2.0 128. Using FishInspector images were assessed for total area (contourDV_regionpropsArea), embryo length (contourDV_regionpropsLengthOfCentralLine), distance between the center of the eyes (YdistanceCenter_eye1DV_eye2DV), and telencephalon distance (YdistanceEdge_eye1DV_eye2DV). Statistical analysis was done in R using the ggsignif package using the Wilcoxon test option 126,129.
Fluorescent images were taken using the Andor Dragonfly High Speed Confocal Platform with the iXon Ultra camera. Human YTHDF2 mRNA (0.5 nL/cell, 100ng/uL human YTHDF2 mRNA) microinjected Tg(HuC-eGFP) strain zebrafish, which harbor a GFP fluorescent pan-neuronal marker were bathed in 0.003% 1-phenyl-2-thiourea (PTU) in 10% Hank’s saline between 20-24pf for 24hrs (Fisher Scientific, 5001443999). At 3dpf zebrafish embryos were embedded in 1% low melt agarose (Thermo Fisher Scientific, BP160-100) and imaged for GFP.
Zebrafish Seizure Analysis
Zebrafish movement was recorded at 30 frames per second for 1 hr with no stimulation using ViewPoint’s ZebraBox Technology per manufacturer recommendations 61. Embryos were maintained at 37C using a polystat water heating system. Movements were detected over 1 second intervals. Distance and time moved were analyzed using a custom R script and high speed events were quantified using a previously published MATLAB script 62.
Zebrafish In Situ Hybridization
Sense and anti-sense probes were designed to ythdf2 and generated via the in vitro transcription kit mMESSAGE mMACHINE™ SP6 Transcription Kit (Thermo Fisher Scientific, AM1340) (F - GAAATTAATACGACTCACTATAGGGctcaacaaggcagtggggat; R - GAAATAATTAACCCTCACTAAAGGGAAgatccagcgcacatcgaaac). NHGRI-1 zebrafish embryos were collected and treated with 0.003% 1-phenyl-2-thiourea (PTU) in 10% Hank’s saline between 20-24pf for 24hrs (Fisher Scientific, 5001443999). At 3dpf PTU treated-embryos were fixed in 4% paraformaldehyde for 12hrs at 4C (Fisher Scientific, AC416785000). Embryos were washed overnight at 70C in a solution containing 100ng of DIG-labeled probes in 250uL Hybridization Mix (Hybridization Mix: 32.5mL 65% deionized formamide, 12.5mL 5X SSC, 0.1% Tween 20, 0.1mL 50ug/mL heparin, 0.25mL 500ug/mL tRNA Type X, 0.92mL 9.2mM citric acid, 3.48mL water) (20X SSC: 87.65g NaCl, 50.26g citric acid trisodium salt dihydrate, to 500mL with nuclease-free water). Embryos were blocked for 4hrs at room temperature (Blocking solution: 0.2mL 2% sheep serum, 2mL 2mg/mL BSA, 7.8mL PBS-Tw, to 10mL with nuclease-free water) (PBS-Tw: 1xPBS containing 0.1% Tween-20) and treated with anti-DIG antibody (1:5000 in Blocking solution) overnight at 4C. Embryos were stained in the dark (staining solution: 33uL NBT and 16.5uL BCIP in 5mL AP buffer) (AP buffer: 100mL 100mM Tris, 20mL 100mM NaCl, 5mL 5mM MgCl2, 2mL 0.1% Tween 20, to 1L with nuclease-free water) for 1hr at room temperature, and destained in methanol. Embryos were visualized via brightfield using a Leica M205 FA microscope.
Data Availability
All data produced in the present work are contained in the manuscript.
ACKNOWLEDGEMENTS
Thank you to the APP research staff, especially Brianna Heath and Melissa Regester, as well as the undergraduate students that provide husbandry and care for our zebrafish. We are grateful to Bruce Appel from the University of Colorado for sharing the HuC-eGFP line. MYD is supported by the NIH (DP2MH119424), SSN is supported by the NIMH Autism Research Training Program T32 (MH073124) through the UC Davis MIND Institute; NAFM is supported by the NIH as a Postbaccalaureate Research Education Program fellow (R25GM116690). A special thank you to the families who have generously shared their genetic data – you are who all of this research is for.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.
- 20.↵
- 21.↵
- 22.
- 23.
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.
- 47.
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.
- 90.
- 91.
- 92.↵
- 93.↵
- 94.
- 95.
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
- 112.↵
- 113.↵
- 114.↵
- 115.↵
- 116.↵
- 117.↵
- 118.↵
- 119.↵
- 120.↵
- 121.↵
- 122.↵
- 123.↵
- 124.↵
- 125.↵
- 126.↵
- 127.↵
- 128.↵
- 129.↵
