Targeted sequencing and integrative analysis to prioritize candidate genes in neurodevelopmental disorders

Neurodevelopmental disorders (NDDs) are a group of diseases characterized by high heterogeneity and frequently co-occurring symptoms. The mutational spectrum in patients with NDDs is largely incomplete. Here, we sequenced 547 genes from 1102 patients with NDDs and validated 1271 potential functional variants, including 108 de novo variants (DNVs) in 78 autosomal genes and seven inherited hemizygous variants in six X chromosomal genes. Notably, 36 of these 78 genes are the first to be reported in Chinese patients with NDDs. By integrating our genetic data with public data, we prioritized 212 NDD candidate genes with FDR < 0.1, including 17 novel genes. The novel candidate genes interacted or were co-expressed with known candidate genes, forming a functional network involved in known pathways. We highlighted MSL2, which carried two de novo protein-truncating variants (p.L192Vfs*3 and p.S486Ifs*11) and was frequently connected with known candidate genes. This study provides the mutational spectrum of NDDs in China and prioritizes 212 NDD candidate genes for further functional validation and genetic counseling. Supplementary Information The online version contains supplementary material available at 10.1007/s12035-021-02377-y.


Introduction
Neurodevelopmental disorders (NDDs) are complex and heterogeneous disorders characterized by impaired motor function, learning, and verbal and non-verbal communication resulting from dysfunction during brain developmental [1][2][3][4]. NDDs span a very wide range of neurological and psychiatric disorders [1,5] and affect more than 3% of children worldwide [6]. Although each diagnosis is distinct in clinical settings, in general, NDDs are characterised by high clinical and genetic heterogeneity. Previous studies have estimated the heritability of NDDs, with the highest being up to ~80% [7][8][9][10]. Multiple phenotype-genotype correlation studies have revealed that patients carrying deleterious variants in the same risk gene have manifested broad clinical phenotypes, including impaired social interaction, developmental delay, regression, and repetitive behaviour [11][12][13]. The high level of heterogeneity and similar pathogenetic mechanisms present challenges for clinical diagnosis and treatment.
Here, we used targeted sequencing to examine 547 genes from 1,102 Chinese patients with NDDs. We identi ed potential functional variants and explored the patterns of DNVs in Chinese patients with NDDs. We then integrated our dataset with public datasets of neurodevelopmental disorders from the Gene4Denovo [29] database to prioritise NDD candidate genes and discover novel candidate genes. Finally, we investigated whether novel candidate genes were functionally associated with known candidate genes using multilevel bioinformatics analysis.
A total of 935 unrelated trios (probands and their unaffected parents) and 167 probands without parents were recruited from China. Genomic DNA (1 μg) extracted from whole blood was sheared and assembled into a DNA library prior to targeted sequencing. The Illumina X10 sequencing system (Illumina, San Diego, CA, USA) was used to generate paired-end raw data. This panel resulted in an average depth of 181.46X in target regions, and 98.82% of target bases were covered with depth ≥ 10X on average. This study was approved by the Institutional Review Board of the State Key Laboratory of Medical Genetics, School of Life Sciences, Central South University, Changsha, Hunan, China. All subjects who participated in this study provided informed consent prior to sample collection.

Variation detection and annotation
Quality control of the sequencing data was performed using Cutadapt [33] and FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to remove adapter and unquali ed sequences, respectively. BWA-MEM [34] was employed to align the clean reads to the human reference genome (hg19).
Samtools [35] utilities were used to mark duplicate reads and generate position-sorted les, while the Genome Analysis Toolkit HaplotypeCaller was used to call variants. Comprehensive annotation of all variants was performed using ANNOVAR [36], including functional implications (gene region, functional effect, mRNA GeneBank accession number, amino acid change, cytoband, etc.), functional predictions for missense variants, and allele frequencies of gnomAD, ExAC, and inhouse data (1,113 WES samples and 2,469 WGS samples). Deleterious missense variants (Dmis) were predicted using ReVe [37] and missense variants with ReVe ≤ 0.7 were excluded. Only protein-truncating variants (PTVs, including stop-gain, stop-loss, frameshift, and splicing) and Dmis with minor-allele frequencies ≤ 0.1% were de ned as potential functional variants and selected for further analysis. Sanger sequencing was used to validate all potential functional variants in our study.

Prioritisation of candidate genes
Using de novo PTVs, together with Dmis and background DNV rates, we employed the TADA classi cation tool [38] to prioritise candidate genes, with gene exhibiting a false discovery rate (FDR) < 0.1 de ned as candidates. First, DNVs from 935 trios, inherited variants from 1,102 probands, and genetic variants from 3,582 inhouse Chinese controls were applied to TADA model and we prioritised 17 candidate genes with FDR values < 0.1. To increase the statistical power, we integrated DNVs from 16,807 probands with different types of NDDs and included 3,391 controls from the Gene4Denovo database (version 1.0) [39] and 208 candidate genes with FDR values < 0.1 for further analysis. All samples integrated in our study have been carefully deduplicated.
To discover novel candidate genes, we collected known candidate genes from previous studies and related databases, and excluded them from the list of candidate genes prioritised by TADA. Known candidate genes were de ned as follows: 1) genes de ned as risk genes in any ten recent publications involving large-scale exome, whole-genome, or targeted sequencing (De Rubeis et al. [17], Sanders et al. [40], Lelieveld et al. [41], McRae et al. [42], Stessman et al. [43], Yuen et al. [44], Nguyen et al. [22], Takata et al. [23], Coe et al. [21], and Satterstrom et al. [45]); 2) genes collected from OMIM that were associated with neurodevelopmental disorders, including intellectual developmental disorder, autism, mental retardation, epileptic encephalopathy, and schizophrenia; 3) genes belonging to the category of syndromic gene or score category of 1 or 2 in SFARI Gene [30]. Genes that met any of the above conditions were classi ed as known candidate genes. In contrast, genes that failed to meet any of the above conditions above and had no obvious genetic evidence in association with NDDs in PubMed were classi ed as novel candidate genes.

Permutation test
Spatial and temporal expression data of the human brain were downloaded from BrainSpan(http://www.brainspan.org/). We calculated the Pearson correlation coe cients between any two genes based on their expression levels. Gene pairs with |R| > 0.8 were regarded as being co-expressed in the human brain, as described in previous study [46]. The permutation test was performed to compare the novel and the known candidate gene sets, in order to evaluate their functional connections. In brief, we compared the number of co-expressed genes within the novel and the known candidate gene sets and their connections with a random gene set of 1,000,000 random iterations.

Functional network analysis
We downloaded PPI data from the STRING database (https://string-db.org/) for analysis and considered gene pairs with a con dence score > 400 to be interacting genes. We then constructed a network based on gene pairs selected from our PPI and co-expression analyses. Novel candidate genes were de ned as 'seed genes' that were directly connected and used to form an interconnected functional network. Known NDD candidate genes directly connected to at least two novel candidate genes were added to the above network. To further investigate the functional pathways of novel candidate genes, we performed GO enrichment analysis using MetaScape (https://metascape.org). Similar pathways were merged into a single cluster. Network gures were drawn using Cytoscape v.3.7.2.
We then employed probability of loss-of-function intolerance (pLI) [47] analysis from ExAC and residual variation intolerance scores (RVIS) [48] to evaluate the functional impact of genes with multiple DNVs. According to the scores of pLI and RVIS, we ranked them and used percentiles as an indicator of gene intolerance. Genes with lowpercentile RVIS or pLI were more likely to be intolerant to genetic variants. As expected, it was shown that genes with multiple DNVs both exhibited signi cantly lower percentile pLI (p = 4.10 × 10 -10 ) and RVIS (p = 1.55 × 10 -10 ) than background genes (Fig. S2a). In particular, 20 of the 21 multiple DNV genes ranked in the top 50% of pLI and RVIS ( Table 1), suggesting that they are less tolerant of damaging variants.

Inherited X-linked hemizygous variants
We identi ed seven inherited X-linked hemizygous variants in six genes (Table S5). Among these seven variants, three variants, including two Dmis (p.R270C on PTCHD1 and p.R197H on SLC16A2) and one splicing variant (c.1030-1G>C on SLC9A7), were recorded in the dbSNP or gnomAD database. The remaining four variants, including three Dmis (p.Q153R on ATP6AP2, p.R287C on ARHGEF9, and p.C517Y on PLXNA3) and one frameshift (p.N136Kfs*30 on SLC16A2) are reported here for the rst time. Compared with the 805 known candidate genes, we note that ve of the above genes are classi ed as known candidate genes and three (PTCHD1, ATP6AP2, and SLC9A7) are reported for the rst time in Chinese patient with NDDs. In addition, we identi ed a novel candidate gene (PLXNA3) carrying a novel hemizygous missense variant (p.C517Y). A previous study has described an ASD patient from Maghreb carrying a rare inherited missense variant (p.D863E) in PLXNA3 [49], suggesting that it may contribute to autism.

Functional characteristics of prioritised candidate genes
We next employed MetaScape [59] to perform functional enrichment analysis. As expected, these candidate genes were signi cantly enriched in NDD-associated pathways, such as synapse organisation, covalent chromatin modi cation, head development, behaviour, and regulation of ion transport [17,[60][61][62][63][64] (Fig. S3). Interestingly, the novel candidate genes were also involved in similar biological pathways. For example, MSL2 was implicated in covalent chromatin modi cation, LRRC4 was involved in chemical synaptic transmission, and PLXNA3 was associated with head development and axonogenesis.
We also performed functional cell-and tissue-speci c enrichment analyses [65,66] to investigate whether candidate genes were associated with speci c tissues or cells. As expected, the 212 NDD candidate genes tended to be enriched in the cerebellum, cortex, and striatum across the early to mid-fetal, mid-late childhood, and young adulthood stages of development (Fig. S4a). In the cell-speci c enrichment analyses, we observed a highly signi cant enrichment in Drd1+ and Drd2+ medium spiny neurons of neostriatum and rods (Fig. S4b), similar to data reported in a previous study [21]. These results suggest that the 212 NDD candidate genes and novel candidate genes are functionally associated with the etiology of NDDs.

Functional network analysis between known and novel candidate genes
To investigate correlations between novel and known candidate genes, we performed a permutation test to estimate the relationship between these genes based on co-expression gene pairs identi ed from the BrainSpan atlas. We found that nine of the 17 novel candidate genes (p = 0.013, permutation test, Fig. S5a) were coexpressed with 149 known candidate genes (p = 0.01, permutation test, Fig. S5b), with 217 connections between them (p = 5.40 × 10 -3 , permutation test, Fig. S5c), suggesting that the novel candidate genes are signi cantly coexpressed with the known candidate genes and are more likely to be related to the pathology of NDDs.
To further investigate the functional relationship between novel and known candidate genes, we constructed a functional network by integrating PPI data from STRING and brain expression data from BrainSpan. Only known candidate genes directly interacting/co-expressed with at least two novel candidate genes were added to the network. The co-expressed/PPI network encompassed 78 genes, including 10 novel candidate genes and 68 known candidate genes (Fig. 2). These genes were enriched in several biological processes known to be related to NDDs, such as covalent chromatin modi cation (GO:0016569, p = 1.55 × 10 -12 , Fisher's exact test), dendritic spine development (GO:0060996, p = 4.74 × 10 -6 , Fisher's exact test), and brain development (GO:0007420, p = 1.92 × 10 -9 , Fisher's exact test). Furthermore, these genes showed a signi cant enrichment of previously reported gene sets: FMRP targets [67] (p = 9.82 × 10 -13 , Fisher's exact test), genes involved in synaptic genes [68] (p = 0.012, Fisher's exact test), and genes essential in mice [69] (p = 1.14 × 10 -18 , Fisher's exact test). It is worth noting that 10 novel candidate genes (UBR3, PLXNA3, ITSN1, POLR3A, MSL2, LRRC4, PSD3, DCX, SPAG9, and SMAD6) were signi cantly involved in functional clusters known to be associated with NDDs. For example, MSL2, UBR3, and PLXNA3 are involved in chromatin organization, while ITSN1, POLR3A, LRRC4, UBR3, and PLXNA3 have been implicated in nervous system development. In addition, we found that 10 novel candidate genes showed coexpression/interaction with more than three known candidate genes (Fig. S6). For example, MSL2, which was the most frequently connected novel gene, was co-expressed/interacted with 35 known candidate genes. DCX, ITSN1, and POLR3A were connected with more than 10 known candidate genes. These results suggest that novel candidate genes are functionally associated with known risk genes and that these 10 novel candidate genes may have a stronger in uence in the etiology of NDDs.

Discussion
Previous studies have provided evidences that targeted sequencing and integrative analyses play a critical role in the discovery of novel candidate genes [21][22][23][24]. Based on targeted sequencing of 547 target genes, we identi ed 108 DNVs, accounting for ~10.91% (102/935) of our Chinese cohort. Among the genes with DNVs in our study, most genes have been reported to be associated with NDDs in previous studies, suggesting that genes identi ed in our Chinese cohort are relevant in other populations. Among the 108 DNVs, 60 DNVs were not recorded in public databases (including dbSNP, gnomAD, Gene4Denovo [39], and PubMed), consistent with the ndings of Georgi et al., who reported that recurrently mutated amino acid sites in genes are rarely detected [70]. These results highlight the role of DNVs in the genetic heterogeneity of NDDs. In addition, we identi ed 21 genes with multiple DNVs present in 51 Chinese patients, and our analyses suggest that these genes display signi cant intolerance to damaging variants. Consistent with previous ASD study [19], SCN2A is the most frequent gene with multiple DNVs in the Chinese population. Compared with the known candidate genes, 20 of 21 genes with multiple DNVs were classi ed as known candidate genes. Notably, this study is the rst to report TCF20 (with its three de novo PTVs) as a candidate NDD gene in a Chinese cohort.
By integrating information regarding DNVs from public datasets with our results, we prioritised 212 candidate genes, con rming that combing public NDD datasets is bene cial for the discovery of candidate genes [22][23][24].
Consistent with our analyses of genes with multiple DNVs, the 212 identi ed candidate genes were more intolerant of damaging variants. In addition, we demonstrated that the 212 candidate genes are closely associated with the etiology of NDDs from the perspective of biological pathway and functional cell-and tissuespeci c enrichment analyses. These results suggest that most of the 212 candidate genes identi ed in our study truly contribute to NDDs and they are worth validating in genetic functional studies or replicating in cohorts.
In this study, we prioritised 17 novel candidate genes and revealed similar functional characteristics between these novel candidate genes and other known candidate genes. Through functional network analysis, we observed that novel candidate genes frequently interacted/were co-expressed with known candidate genes, and genes in the network were enriched in NDD-associated clusters, as described in previous studies [4,17,[71][72][73].
Interestingly, 10 novel candidate genes were closely connected with known candidate genes and were involved in NDD-associated clusters, suggesting that these novel candidate genes are more likely to be associated with NDDs. For example, SMAD6 (with an FDR value < 0.001 in the TADA analysis)belongs to an essential gene set in our study and was connected with ve known candidate genes. Timberlake et al. reported a role for SMAD6 in craniosynostosis and a high risk of developmental delay [74]. Furthermore, in their functional study, Xie et al.
showed that Smad6 inhibits the Wnt/beta-catenin pathway and plays a role in neuronal differentiation in the intermediate zone of the chick developing dorsal spinal cord [75]. Our analyses revealed that SPAG9 carried a de novo PTV in our cohort and it is reportedly overexpressed in human astrocytoma which arises from neural progenitor cells in the central nervous system [76]. We also found that SPAG9 was implicated in chromatin organisation and co-expressed/interacted with ve known candidate genes, further supporting a role for this gene in NDDs. Among the novel candidate genes, MSL2 was a particularly interesting gene that carried two de novo PTVs in our Chinese cohort. Except for the multiple DNVs in MSL2, it was the gene that most frequently coexpressed/interacted with known candidate genes and was involved in the functional cluster of chromatin organisation. Iossifov et al. reported a de novo MSL2 PTV in a patient with autism [71], further suggesting that MSL2 may be strongly linked to NDDs. We expect to have a large cohort study or functional experiments to validate this possibility in future studies.
Despite our best efforts to comprehensively integrate public data with our cohort data to discover novel risk genes, our study was still bound by some limitations. First, while we discovered that novel candidate genes identi ed in our study are functionally associated with known candidate genes at multiple levels, these novel candidate genes lack strong evidences that support their roles in NDDs. We expect that these novel candidate genes could be validated in larger cohorts or through functional genetic experiments. Second, the sample size of our Chinese cohort is not large enough, and future studies need to recruit more volunteers to improve the statistical power and allow for comparisons with variant patterns of other populations. Third, since NDDs are characterised by high clinical and genetic heterogeneity, and multiple risk factors contribute to the etiology of these disorders, candidate gene analyses can only partially elucidate the processes underlying NDDs. Further studies are needed to integrate the impact of distinct in uences such as epigenetic, environmental, and genetic factors.
In summary, our study demonstrates that integrating DNVs from multiple NDD-related studies can help in the identi cation of risk genes. We highlight that the pattern of DNVs in Chinese cohorts is relevant to other populations. Furthermore, we provide evidence at multiple levels that novel candidate genes are functionally associated with known candidate genes. Finally, our study describes new high-con dence risk genes that should aid the study of the NDD etiology and we expect that these genes will be worth analysing in large cohorts or being