Rare genetic variants confer a high risk of ADHD and implicate neuronal biology

The iPSYCH study was approved by the Scientific Ethics Committee in the Central Denmark Region (case number 1-10-72-287-12) and the Danish Data Protection Agency. In accordance with Danish legislation, the above-mentioned ethics committee waived the need for specific informed consent in biomedical research based on existing biobanks. iPSYCH was approved by the ethics committee in 2012, with subsequent amendments in 2013, 2015 and 2018. More details can be found at https://ipsych.dk/en/data-security/health-research-and-ethical-approval. New Danish legislation (effective from January 2024) introduces the possibility for participants to opt out of studies that are exempt from active informed consent. After consulting with the Ethics Committees and patient organizations, iPSYCH contacted all participants (around 140,000) in the iPSYCH cohort in June 2025 and offered the possibility of opting out of new genetic studies initiated henceforth. Overall, 1.8% of the iPSYCH participants chose to opt out, and their data will be deleted from the active research database. Data included in finalized and ongoing studies will not be removed.

The clinical data were approved by the ethics committee at the University of Würzburg in Germany. In the Netherlands, they were approved by the regional ethics committee (Commissie Mensgebonden Onderzoek: CMO Regio Arnhem—Nijmegen; protocol III.04.0403 and 2014/290; ABR: NL47721.091.14) and the Institutional Review Board of the Radboud University Medical Center. Participants were included at the Department of Psychiatry at the Radboud University Nijmegen Medical Centre. All participants in the German and Dutch samples provided signed informed consent in accordance with the Declaration of Helsinki.

Samples

iPSYCH

The individuals selected for exome sequencing were part of the iPSYCH cohort, which has been described in detail elsewhere^15,16. In short, the study base includes all singleton births to mothers who were living in Denmark between 1 May 1981 and 31 December 2008, where the child was alive and resided in Denmark at their one-year birthday (n = 1,657,449). All individuals diagnosed with major psychiatric disorders by the end of 2016 according to the ICD10 criteria (ADHD (1.8% in the study base), autism, bipolar disorder, schizophrenia, major depressive disorder or post-partum depression (n = 93,608)) were identified in the study base using information in the Danish Psychiatric Central Research Registry⁴⁹ (and the Danish Patient Registry⁵⁰ for some disorders). In addition, 50,000 randomly selected population-based controls from the study base were selected. Subsequently, biological material for genotyping was obtained from the Danish Neonatal Screening Biobank (DNSB)⁵¹. The DNSB has stored residual biological material from screening of newborns for rare metabolic disorders since May 1981, and includes material from practically all births in Denmark since then. A subsample of 34,544 individuals was selected for whole-exome sequencing, and from these, we included individuals with an ICD10 diagnosis of ADHD (F90) in the Danish Psychiatric Central Research Registry⁴⁹ and the Danish Patient Registry⁵⁰ given before or during 2016; individuals with no diagnosis of the major psychiatric disorders (autism (ICD10: F84.0, F84.1, F84.5, F84.8 and F84.9), bipolar disorder (ICD10: F30–F31), schizophrenia (ICD10: F20) or major depressive disorder (ICD10: F32–F33)) were included as controls. The samples were included in iPSYCH in 2012–2016 and the sequencing was performed in 2012–2018.

For the comorbidity analyses, we identified individuals with the following ICD10 diagnosis codes in the Danish Psychiatric Central Research Registry: ID (ICD10: F70, F71, F72, F73, F78 and F79); ASD (as above); schizophrenia (as above); DBDs (including conduct disorder and oppositional defiant disorder (ICD10: F91 and F90.1); SUDs (ICD10: F10.1–9, F11.1–9, F12.1–9, F13.1–9, F14.1–9, F15.1–9, F16.1–9, F17.1–9, F18.1–9 and F19.1–9) and multi-comorbidities, which, besides the comorbidities already listed, included comorbid anxiety (ICD10: F40.0–F40.2, F41.0–F41.1, F42 and F43.0–F43.1), tic disorder (ICD10: F95), bipolar disorder (ICD10: F30–F31), major depressive disorder (ICD10: F32–F33), anorexia nervosa (ICD10: F50.0), DDs (ICD10: F80–F83) and antisocial personality disorder (ICD10: F60.2).

Clinical samples

The clinical samples consisted of adults (over 18 years old) with persistent ADHD diagnosed according to the Diagnostic and Statistical Manual of Mental Disorders IV (DSM-IV) criteria. None of the individuals was diagnosed with ID. They were recruited as part of the International Multicenter Persistent ADHD Collaboration (IMpACT) at two sites: Radboud University Medical Center, the Netherlands, and University Hospital Würzburg, Germany. In analyses of these clinical samples, we used control samples from individuals recruited together with the clinical cases at the IMpACT site at Radboud University Medical Center (ADHD-screened controls), and from 1,766 German control individuals who were recruited from the German MI Family Study⁵² and the Angio-Lub study; the latter samples were whole-exome-sequenced by the MIGen Exome Sequencing Consortium: Lubeck Heart Study (dbGaP accession number phs000990/DS-CVD, https://dbgap.ncbi.nlm.nih.gov/beta/study/phs000990.v1.p1/). The German controls consisted of 870 individuals with cardiovascular disease and 896 without. The total set of samples described in this section is referred to as the ‘clinical sample’ in the remainder of this manuscript.

GnomAD

All references to gnomAD refer to release 2.1.1 exomes from a subset of gnomAD consisting of individuals with non-Finnish European ancestry, and no diagnosis of psychiatric and neurological disorders (n = 44,779) (see ‘Data availability’).

Exome sequencing and quality control

iPSYCH

In this study, we applied the same methods for quality control (QC) described in our previous study¹⁴ to an updated dataset including new exome-sequenced individuals. To recap, the sequencing of 34,544 individuals from the iPSYCH cohort was performed using the Illumina Nextera Capture kit and the Illumina HiSeq sequencer. The sequencing process was performed in waves, including a pilot wave and three more substantial production waves. After sequencing, raw data were processed using the Genome Analysis Toolkit (GATK) v.3.4 to generate a variant call format (VCF) v.4.1 file. As per Danish regulations, American College of Medical Genetics (ACMG) v3.2 genes⁵³ were removed from the dataset. Next, we did thorough, multiple-round quality checks on the samples and the genetic variations using Hail 0.1. Samples were removed if they lacked complete phenotype information; if the imputed sex was inconsistent with the reported sex in the registries; if duplicates or genetic outliers were identified by principal component analysis (PCA); if they had an estimated level of contamination greater than 5%; or if they had an estimated level of chimeric reads higher than 5%. In addition, one of each pair of related samples was removed if the pairwise pi-hat value was 0.2 or higher.

Genotypes were removed if they did not pass GATK variant quality score recalibration (VQSR) or had a read depth lower than 10 or higher than 1,000. QC was done in Hail 0.1. Homozygous alleles were removed if they had reference calls with a genotype quality lower than 25, homozygous alternate alleles with PL(HomRef) (that is, the phred-scaled likelihood of being homozygous reference) < 25 or less than 90% of reads supporting the alternate allele. Heterozygous alleles were removed if they had PL(HomRef) < 25 or less than 25% of reads supporting the alternate allele, less than 90% informative reads (that is, number of reads supporting the reference allele plus number of reads supporting the alternate allele < 90% of the read depth) or a probability of the allele balance (calculated from a binomial distribution centred on 0.5) less than 1 × 10⁻⁹. After these filters were applied, variants with a call rate lower than 90% were removed, then samples with a call rate lower than 95%, and then variants with a call rate lower than 95% were removed. After QC, 28,448 individuals and 1,362,971 variants remained for further analysis.

Subsequently, we selected for this study the individuals diagnosed with ADHD and controls as described above, resulting in 8,895 cases and 9,001 controls (see Supplementary Table 23 for sample of individuals with other ancestries not included in this study). Finally, we defined rare variants (n = 565,053) as those with an allele count no higher than five across our dataset (n = 17,896) and the exome subset of gnomAD used in this study (n = 44,779).

Clinical samples

Biological samples collected at the Dutch and German IMpACT sites were sequenced at BGI, Shenzhen, China. The coding regions of the DNA were targeted using BGI´s exome-capture kit (developed by Beijing Genomics Institute targeting 58.8 Mb) and paired-end sequenced on the Illumina HiSeq 2000 platform, with an average sequencing depth of 50×. The exome-sequencing data obtained from dbGaP were generated at the Broad Institute of Harvard and MIT using Illumina’s ICE Capture reagent kit, and sequencing was performed on an Illumina HiSeq 2000 or 2500. Bam files for all samples (clinical cases and dbGaP controls) were reprocessed together, adhering to the same QC approach described above for the iPSYCH samples. After this QC process, 2,816 samples, consisting of 1,078 individuals with ADHD and 1,738 controls, remained for further analysis. We were not able to obtain a group of completely homogenous individuals after the exclusion of genetic outliers based on PCA (Supplementary Fig. 11), and thus, to minimize the effect of population stratification in subsequent analyses, we aimed at getting as close to de novo mutations as possible; we therefore included only ultra-rare variants, defined as singletons not in the non-psychiatric subsample of individuals in the ExAC¹⁸ database. We could not filter based on presence in gnomAD because the dbGaP data (that is, the controls) are included in this database.

Effects of variant categories on ADHD

Quality controlled variants were functionally annotated using SnpEff v.4.3^54,55, and SnpSift⁵⁴ was used to annotate information derived from dbNSFP⁵⁶. If a variant had more than one annotation, only the most severe annotation was considered. PTVs were defined by being annotated as frameshift, splice-site or stop-gained, and predicted with a loss of function (LOF) flag by SnpEff. For missense variants, we annotated the potential effect of the variant on protein function using the MPC score²⁰.

In iPSYCH data, the burden of different rare-variant categories (rPTVs; rSevereDMVs, rModerateDMVs and rSYNs) in ADHD cases compared with controls was tested using logistic regression. Following the same approach, the load of rare variants across variant categories in female individuals with ADHD (n = 2,265) was compared with male individuals with ADHD (n = 6,630), both with and without individuals with ID (1,827 female individuals without ID; 5,414 male individuals without ID). The burdens in all genes and genes stratified by their pLI score were evaluated. Covariates included in analyses of iPSYCH samples were: birth year, sex, the first ten principal components (PCs) from PCA (performed after excluding non-European samples), total number of variants, number of rare synonymous variants, percentage of exome target covered at a read depth of at least 20, mean read depth at sites within the exome target passing VQSR, and sequencing wave (one-hot encoded).

The burden of different types of variant categories in gene sets was also tested for the clinical samples following the approach described above but restricting to ultra-rare variants and including the following covariates: sex, total number of variants, number of ultra-rare synonymous variants, number of all variants and PC1–PC10 from PCA (performed after excluding individuals with non-European genetic ancestry).

Gene-based burden tests

To increase power for gene discovery, we combined iPSYCH data with a subset of 44,779 individuals from gnomAD¹⁷ (as defined above). Only variants in high-confidence regions for the two datasets were included, defined as regions in which at least 80% of the samples in both datasets had at least 10× sequencing coverage (based on analysis of bam files for the Danish samples and on coverage summary tables for gnomAD). To avoid biases caused by variations in call rates between cases (entirely iPSYCH) and controls (83.3% gnomAD), all autosomal genes with higher rates of rare synonymous variants in cases than controls were excluded (15,603 out of 18,866 autosomal genes remained; 3,263 excluded).

Variants were grouped depending on their impact on ADHD: class I variants include rPTVs and rare missense variants with MPC > 3 (rSevereDMVs); class II variants include missense variants with 2 ≤ MPC ≤ 3 (rModerateDMVs). Gene-based burden analysis was performed with a two-tailed Fisher’s exact test. Only genes with an increased load of class I (3,698 genes) and/or class II (1,026 genes) variants in cases compared with controls were considered. We did not consider genes with higher numbers of class I and II variants in controls, because such observations could be caused by a generally higher number of variants in controls, owing to the previous filtering step in which only genes with a higher number of rare synonymous variants in controls than in cases were retained. A small number of genes (347) had an increased load of both class I and class II variants in cases, and for these genes we also estimated the combined impact of class I and class II variants in a meta-analysis using the weighted z-score method. Weights were the ratio of the standardized effect sizes observed for the classes in enrichment analysis of constrained genes. For the 347 genes, we used the minimum P value across the three analyses done for these genes. In total, we tested 15,603 genes (347 genes were tested three times) resulting in 16,297 tests. We considered a gene significantly associated with ADHD if the P value was lower than 0.05/16,297 = 3.07 × 10⁻⁶.

As described above, genes with an increased rate of rare synonymous variants in individuals with ADHD relative to controls were excluded to avoid biases when combining iPSYCH and gnomAD data. When inspecting the QQ-plot of rare synonymous variants (Supplementary Fig. 12), the plot revealed the anticipated pattern of inflation that corresponds to our gene selection strategy.

In the clinical sample, we performed a gene-based burden test of class I variants using emmaxCMC^57,58 implemented in EPACTS v.3.3.0 (https://genome.sph.umich.edu/wiki/EPACTS). This method allows for the incorporation of a kinship matrix to account for potential remaining population stratification among individuals. We generated the kinship matrix using ‘epacts make-kin –min-maf 0.01 –remove-complex’. In addition, we included sex, the number of ultra-rare synonymous variants, the number of all variants and PC1–PC10 from PCA (performed after excluding non-European samples) as covariates.

Gene-set analyses in the clinical sample

We tested for increased load of ultra-rare class I and class II variants in ADHD compared with controls in the clinical sample. We used different P-value thresholds from the gene-based burden test in iPSYCH + gnomAD samples to define three ADHD rare-variant risk gene sets: P < 1 × 10⁻³ (20 genes), P < 5 × 10⁻² (316 genes) and P < 1 × 10⁻¹ (583 genes) (gene sets are listed in Supplementary Table 20). The control group included 870 with cardiovascular disease and 896 without, so we also tested for the potential difference in ultra-rare class I variants across the two groups in the three ADHD gene sets. Only class I variants were tested, because these were the only type of variants with a tendency to be overrepresented in clinical ADHD cases versus controls (Supplementary Information, Supplementary Fig. 13 and Supplementary Table 24). The gene sets were tested using logistic regression with sex, number of ultra-rare synonymous variants, number of all variants, and PC1–PC10 from PCA (performed after excluding non-European samples) as covariates.

X-chromosome analyses

Given the distinct regions of the X chromosome, we applied region-specific thresholds to define rare variants: in the pseudoautosomal regions (PARs), rare variants were defined as those with an allele count of five or less following what was done for the autosomes. In non-pseudoautosomal regions (nonPARs), the rare-variant threshold was adjusted to an allele count of three of less to account for the hemizygosity in male individuals. This adjustment was derived using a scaling factor based on the relative contributions of male individuals and female individuals in iPSYCH (6,002 female and 11,894 male individuals) and subset of gnomAD used in this study (19,916 female and 24,863 male individuals):

Variants were categorized into class I, class II and synonymous variants, using the same criteria applied for autosomal variants. In the burden test of male individuals, the number of alleles for nonPAR variants was counted as two to account for hemizygosity.

Logistic regression was performed on three gene sets stratified by their pLI score, comparing the following groups of individuals: ADHD versus controls; ADHD without ID versus controls; ADHD with ID versus controls; and ADHD with ID versus ADHD without ID. Covariates included birth year, sex, the first ten PCs from PCA (non-European samples excluded), total variant count, rare synonymous variant count, exome target coverage (percentage ≥ 20×), mean read depth at target sites and sequencing batch (one-hot encoded). In addition, sex-stratified analyses were performed using the same approach.

For gene discovery on the X chromosome, we applied the same approach as for autosomes when combining iPSYCH data with gnomAD. We excluded genes with higher numbers of synonymous variants in individuals with ADHD than in controls and retained genes with a higher rate of class I variants in cases than in controls (78 genes included) or a higher number of class II variants (44 genes). We performed gene-based burden analysis separately for male individuals and female individuals, in which we, for each gene, compared the number of individuals with ADHD carrying at least one class I or class II variant to the number found for control individuals, using a two-tailed Fisher’s exact test. Notably, six genes had both increased class I and II variants in cases compared with controls; for these genes, a meta-analysis was used to combine the effect of class I and class II variants using the weighted z-score method. The weight was determined as the ratio of the standardized effect sizes observed for class I and class II variants in the enrichment analysis of constrained genes on chromosome X (weight = 1.03/1.86 = 0.55). In total, we performed 78 × 3 tests for class I variants (female, male and combined), 44 × 3 tests for class II variants and 6 tests for combining class I and II, resulting in a total of 372 tests. The Bonferroni correction threshold for statistical significance was set at P < 1.34 × 10⁻⁴.

Rare burden heritability

We used BHR²⁶ to estimate the burden heritability and the contribution of a gene set to the burden heritability (namely, the burden heritability enrichment of a gene set). We used BHR to estimate the heritability of ADHD explained by the load of rare class I, class II and synonymous variants.

Variant-level summary statistics associated with ADHD, including allele count and allele frequency from the iPSYCH exome data, were used as input for BHR. The method regresses gene burden test statistics (based on the variant category being evaluated) against burden scores that correspond to the combined allele frequency. The slope of the regression represents the burden heritability and confounding factors such as population stratification, are controlled through the intercept.

PPI-network analyses

The ANK2 PPI networks were derived from published IP–MS datasets included in Table S2 of a previous study²⁷. All significant proteins with log₂-transformed fold change (FC) > 0 and FDR ≤ 0.1 in the ANK2_WH, ANK2_CNCR1 and ANK2_CNCR2 datasets were defined as the ANK2 interactors in ExNs; the significant proteins in the ANK2_WT dataset were defined as the ANK2 interactors in NPCs.

The MAP1A and ANO8 PPI networks in NPCs and ExNs were derived from IP–MS experiments performed in this study. We evaluated the expression of cell-type marker genes using single-nucleus RNA-seq (snRNA-seq) and performed immunofluorescence staining on NPCs and fully differentiated ExNs to confirm their identities as immature neural progenitors and upper-layer prefrontal cortex neurons, respectively (results shown in Supplementary Fig. 14). Details on cell culture and differentiation, snRNA-seq, immunofluorescence, protein extraction, immunoprecipitation, immunoblotting, mass spectrometry and IP–MS data analysis can be found in the Supplementary Information. Consistent with the ANK2 networks, we defined significant proteins with log₂-transformed FC > 0 and FDR ≤ 0.1 in each IP–MS experiment to be the interactors of the index protein. The resulting IP–MS datasets are provided in Supplementary PPI Tables 1–4.

We parsed the interactors identified across all IP–MS datasets into 12 PPI networks grouped by index proteins and cell-type specificity (Supplementary PPI Table 5). For each index protein, we merged all interactors identified in the same cell type into an ‘NPC’ or ‘ExN’ network, or in either cell type, into a ‘Union’ network. We also merged interactors for all three index proteins into combined NPC, ExN and Union networks accordingly.

We also annotated unique interactors identified across all IP–MS datasets with the following information (Supplementary PPI Table 6): (1) name and number of associated baits (index proteins) in NPCs, ExNs or either cell type; (2) whether the interactor had been implicated in genetic association studies of ADHD (P < 0.001 in this study, or by common-variant risk genes listed in Supplementary Table 7 of Demontis et al.¹³), autism, DDs, neurodevelopmental disorders (Supplementary Table 11 of Fu et al.²⁸) or schizophrenia (Supplementary Table 5 of Singh et al.¹⁹ at various significance thresholds).

For subsequent network enrichment analyses, we compared the network genes against a background set of protein-coding genes expressed in neurons (hereafter, ‘neuronal background’). To define the neuronal background, we re-analysed RNA-seq data derived from the same ExN cellular model used in this study (day-21 and day-51 data from the Gene Expression Omnibus (GEO) GSE178896 dataset)²⁷. We first performed transcript quantification from FASTQ files using Salmon (v.1.10.2)⁵⁹ and GENCODE (v.43)⁶⁰ reference files. We summarized the quantification results to gene-level counts using tximport (v.1.26.1), then removed non-protein-coding genes and low-count genes using the filterByExpr function in edgeR (v.3.40.2). This resulted in a list of 13,018 neuronal background genes (Supplementary PPI Table 7) to be used in downstream analyses.

To perform rare-variant enrichment analysis for the PPI networks, gene-based association scores were obtained from exome-sequencing studies of autism and DDs (FDR_TADA_ASD and FDR_TADA_DD columns in Supplementary Table 11 of Fu et al.²⁸) and schizophrenia (‘P meta’ column in Supplementary Table 5 of Singh et al.¹⁹). Loss-of-function constraint scores (pLI scores) were obtained from gnomAD v.2.1.1¹⁷. For each phenotype and each PPI network, we performed a one-tailed Kolmogorov–Smirnov test to assess whether the network genes had more significant scores than the neuronal background genes. Because the PPI networks were significantly enriched for loss-of-function constrained genes, we also repeated the analysis for constrained (pLI ≥ 0.9) and non-constrained (pLI < 0.9) genes separately. That is, we used one-tailed Kolmogorov–Smirnov tests to compare the network genes with pLI ≥ 0.9 to other neuronal background genes with pLI ≥ 0.9, and vice versa for the network genes with pLI < 0.9.

A description of the common-variant enrichment analysis can be found in the Supplementary Information.

Enrichment analyses

We tested for enrichment of top associated rare-variant risk genes (20 genes with P < 1×10⁻³ in ADHD iPSYCH cases versus iPSYCH + gnomAD controls) among genes in the following gene sets: (1) gene sets related to gene ontology analysed using the Gene Ontology (GO) knowledgebase^61,62 (GO biological processes v.2022-07-01, 9,290 gene sets; GO cellular component ontology v.2022-07-01, 1,581 gene sets; GO molecular function ontology v.2022-07-01, 2,997 gene sets); (2) genes related to synapse function using synaptic annotations based on published, expert-curated evidence for 1,602 genes in SynGO²⁹ (v.20231201 release 1.2); the list of brain-expressed genes provided by SynGO was used as background; (3) gene sets related to biological pathways (BioCarta 2016, 237 gene sets; KEGG 2021, 320 gene sets; Reactome 2022, 1,818 gene sets; canonical pathway gene sets derived from the WikiPathways pathway database 2021, 622 gene sets)²⁹; and (4) gene sets related to diseases (PheWeb v.2019, 1,161 gene sets; PhenGenI Association v.2021, 950 gene sets; GWAS catalogue v.2021, 17,37 gene sets; DisGeNET v.6.0, 9,828 gene sets; OMIM Disease, 90 gene sets). Enrichment analyses of the latter two (that is, pathway- and disease-related gene sets) were performed using Enrichr⁶³. All of the enrichment analyses were done using a one-sided Fisher´s exact test and the within-database correction for multiple testing was done using the Benjamini–Hochberg method. A gene set was considered significant if the within-database corrected P value was less than P = 0.0038 (0.05 divided by the number of databases [0.05/13 = 0.0038]).

In addition, the genes encoding proteins in the union PPI networks (NPCs + ExNs) for each of the three index proteins (MAP1A, ANO8 and ANK2) were tested for enrichment among gene sets related to (1) gene ontology, (2) synapse function and (3) biological pathways, as described above.

Expression of risk genes in the brain

Expression of the top 20 rare-variant risk genes (P < 0.001) across neocortex developmental stages was evaluated using bulk RNA-seq data v.10 (gene-level RPKMs) obtained from BrainSpan (www.brainspan.org). The samples represent an age span from post-conceptual week 8 to 40 years of age, and were grouped into brain developmental stages as defined previously⁶⁴. Following a previous study⁴⁵, we analysed neocortical regions. Samples with poor quality (RIN ≤ 7) were removed. Genes were defined as expressed if the RPKM was at least one in at least 80% of the samples for at least one neocortical region in one major temporal epoch.

After filtering, the BrainSpan dataset contained expression data from 324 samples (information about the number of samples analysed for each developmental stage can be found in Supplementary Table 12). We focused on the 20 ADHD risk genes, of which 17 had expression data in BrainSpan. Gene expression was natural logarithm-transformed (log(RPKM + 1)) and a two-sided paired t-test was used to test for differential expression of the set of the 17 ADHD risk genes against a background gene set (22,402 genes or transcripts) for all developmental stages merged and across developmental stages (P = 4.17×10⁻³ was considered significant correcting for 12 brain developmental stages). A two-sided paired t-test was also performed to determine differential expression of the 17 genes prenatally and postnatally. Likewise, each of the three exome-wide-significant genes (MAP1A, ANK2 and ANO8) were tested for differences in pre- and postnatal expression. Because the sample sizes were small for some developmental stages, we tested to see whether the data were normally distributed. If this assumption was violated, a non-parametric Wilcoxon signed-rank test was performed instead.

Linking risk genes to cell types

We used scDRS³⁰ to link rare-variant risk genes to cell types. First, we used a hypothesis-based approach that focused on midbrain neurons, motivated by findings linking common variants[13] to midbrain dopaminergic neurons. We used four datasets from two published snRNA-seq and scRNA-seq studies (see ‘Data availability’): one scRNA-seq dataset from a study (study I) of gene expression in prenatal human brain cells³¹ derived from ten human embryos 6–11 weeks old (1,695 cells analysed); and three snRNA-seq datasets from a study (study II) of gene expression in developing midbrain neuronal cell types³². These data were generated from 215 pluripotent stem cell (iPS cell) lines, each derived from a single healthy donor (88 male and 127 female individuals), for differentiation towards midbrain neuronal cell types. The three datasets represent cells developed for 11 days (253,381 cells analysed), 30 days (250,923 cells analysed) and 52 days (303,856 cells analysed). We used processed scRNA-seq and snRNA-seq count data and cell-type annotations generated as described in the two papers^31,32; for study I, the data were downloaded as a collapsible expression format (CEF) file and converted into a hierarchical data format version 5 annotated data (h5ad) file format using the Python package AnnData (see ‘Data availability’); for study II, h5ad file formats were available. For study II, the 11-day and 30-day datasets (downloaded) were used without any modifications, whereas the 52-day dataset was filtered to remove cells that were treated with rotenone (following what was done in the published study³²). In addition, the raw read counts from study I (ref. ³¹) and the 52-day data from study II (ref. ³²) were filtered so that only genes represented in a minimum of 30 cells were included.

Furthermore, we used a hypothesis-free approach, in which we evaluated the scDRS in 382 cell-type clusters (2,480,956 cells analysed) representing a spectrum of neurons from the entire human brain (study III; ref. ³³). This dataset consists of snRNA-seq data of neurons desiccated from four adult post-mortem human brains (three male individuals and one female individual) from 105 locations across the forebrain (cerebral cortex, hippocampus, cerebral nuclei, hypothalamus and thalamus), midbrain and hindbrain (pons, medulla and cerebellum). We used a h5ad file with processed data (see ‘Data availability’) and cell cluster definitions reported previously³³.

For the dataset from study I (ref. ³¹) and the 52-day dataset from study II (ref. ³²), there were no generated UMAPs for visualizing cell-type clusters. To generate UMAPs for these data, we did the following, using Python and the single-cell analysis libraries Scanpy and AnnData: for each of the four scRNA-seq datasets, the raw counts were filtered to contain genes represented in a minimum of 30 cells. Afterwards, the raw counts were normalized and log-transformed, and highly variable genes were identified using ‘highly_variable_genes’ in Seurat⁶⁵ implemented in Scanpy with default settings. The expression data for the highly variable genes were scaled and used as input in PCA generated using the singular value decomposition (SVD) solver method in the ‘arpack’ algorithm implemented in Scanpy. PCs from the PCAs were then used to compute a neighbourhood graph generated using the Scanpy function ‘pp.neighbors’, which afterwards was embedded and visualized as a UMAP constructed with the Scanpy function tl.umap. For study I, the neighbourhood graph was constructed with the first 20 PCs and ‘n_neighbors’ set to 100. For study II, the neighbourhood graph was constructed with the first 40 PCs and ‘n_neighbors’ set to 30 for each of the 3 datasets. The parameters were set based on visual inspections.

We used the P values from the top 100 most associated genes from the gene-discovery analysis of iPSYCH + gnomAD samples as input for scDRS analysis, which was done separately for the above-described data from studies I, II and III. The number of genes was set to 100, which is the recommended minimum number of genes for the method. The method computes, on the basis of the expression of these 100 most significant genes, one raw disease score per cell. In addition, 1,000 raw control scores were computed for each cell. The control scores were generated using Monte Carlo sampling, in which 1,000 gene sets were produced, each consisting of 100 genes with similar gene size, mean gene expression and expression variance to that observed for the disease gene set. The computation of scores was corrected for relevant covariates (see below) and the number of genes in each of the cells. The normalized disease score was compared with the empirical normalized control score distribution to estimate a P value that quantifies the association between the disease genes and their expression in individual cells. These scores were used in the downstream scDRS cell-type-level analysis to test the association between the scores and predefined cell types using a t-test and the top 5% quantile of the disease score of cells from the given cell type. In addition, a heterogeneity test using the Geary’s C statistic⁶⁶ was performed to determine heterogeneity in disease score within cell types. Both the compute-score and the perform-downstream functions were run with default settings on raw counts.

We included donor ID as a covariate in analyses of studies II and III. Information on donor ID was not available for study I, and we were therefore not able to correct for the potential effects of factors captured by the donor ID covariate; that is, differences in biological variance between individuals, potential technical variation linked to donor samples or unequal donor representation among cells. The day of cell collection (collected from week 6 to week 11) was also used as a covariate in analyses of data from study I. Sex was used as a covariate in the analyses of data from study II. Information about the sex of the donors was obtained from the Human Induced Pluripotent Stem Cell Initiative (HipSci) data browser (see ‘Data availability’) and merged with the downloaded data. Results were considered significant after within-study Bonferroni correction, correcting for the number of cell types analysed.

Gene-set analyses across comorbidities

The loads of class I and class II variants in individuals with ADHD comorbid with other disorders (ID, autism, schizophrenia, DBDs, SUD and multi-comorbidities; see diagnosis codes above), compared with ADHD without comorbidities, were evaluated in iPSYCH samples for sets of autosomal genes grouped by their pLI score: pLI ≥ 0.9 (2,811 genes), 0.5 < pLI < 0.9 (1,332 genes), pLI ≤ 0.5 (14,267 genes), and for 7 gene sets related to other psychiatric disorders and ID: (1) ‘SCZ_Pval2.14e-6’ includes the 10 genes significantly associated with schizophrenia at P < 2.14 × 10⁻⁶, identified from exome-sequencing data¹⁹; (2) ‘SCZ_Qval0.05’ includes a broader set of 32 genes with a nominal association with schizophrenia at q < 0.05, identified from exome-sequencing data¹⁹; (3) ‘ASD_Pval2.5e-6’ includes 60 genes significantly associated with ASD at P < 2.5 × 10⁻⁶, identified from sequencing data⁶⁷; (4) ‘ASD_FDR0.001’ includes 72 genes significantly associated with ASD at FDR < 0.001, identified from exome-sequencing data²⁸; (5) ‘ASD_FDR0.05’ includes 183 genes with a nominal association with ASD at FDR < 0.05 based on exome-sequencing data²⁸; (6) ‘NDD_FDR0.001’ includes 373 genes significantly associated with neurodevelopmental disorder (NDD) at FDR < 0.001, identified from exome-sequencing data²⁸; and (7) ‘NDD_FDR0.05’ includes 662 genes with a nominal association with NDD at FDR < 0.05, identified from exome-sequencing data²⁸ (gene sets are listed in Supplementary Table 20). The analyses were done using logistic regression and the same covariates as were used in analyses of the effects of variant categories on ADHD in iPSYCH samples (described above). The load of synonymous variants was also evaluated as a sanity check, because we expected no different load for this variant category. We considered tests with P < 0.05/7 = 7.14×10⁻³ as significant (seven gene sets tested).

Effects of rare variants on education and SES

A copy of the whole-exome-sequencing data (after QC and functional annotation) on iPSYCH ADHD cases was transferred to Statistics Denmark (no controls due to restriction on file size), to link rare variants to variables only available at the secured servers at Statistics Denmark. In these analyses, rare variants were defined as singletons in iPSYCH ADHD cases. Individuals with ADHD with at least one rPTV (nonsense, frameshift and essential splice-site variants) in constrained genes (pLI ≥ 0.9) were compared with individuals with ADHD without rPTVs in constrained genes; rPTVs were used as exposure, and high versus low SES or education were used as outcomes in logistic regression analysis using the glm() function in R and for two-tailed Fisher’s exact test using the fisher.test() function in R⁶⁸. We corrected for gender, birth year, first ten PCs generated from ancestry PCA, number of rare synonymous variants, percentage of target with coverage greater than 20×, mean read depth at sites within the exome target passing VQSR, total number of variants, and sequencing wave. As control tests, we analysed rare synonymous variants in constrained genes.

SES information was obtained from the Income Statistics Register⁶⁹; data were available for 5,297 individuals with ADHD who were over 16 years old. The ‘low SES’ group was defined as individuals receiving social security payment, receiving early retirement benefit and/or having been unemployed for more than six months (n = 3,110); the ‘high SES’ group consisted of the remaining individuals (n = 3,223). Education level was obtained from the Danish Population Education Register⁷⁰. Low education was defined as finishing only primary school, which is nine years of schooling in Denmark (n = 6,488); high education was defined as an education beyond primary school (n = 1,463).

We used two-tailed Fisher´s exact test to estimate the risk of low SES or low education among individuals with ‘ADHD with one or more rPTVs on constrained genes’ or individuals with ‘ADHD without rPTVs in constrained genes’ against 21,413 population-based controls (over 16 years old).

We tested for associations of ultra-rare deleterious variants with measures of IQ in adults with ADHD in the German clinical sample (n = 962 individuals). Intellectual function was assessed with the Mehrfachwahl Wortschatz Intelligenztest (MWT-B)⁷¹ test, and none of the recruited individuals had an IQ lower than 80. We used linear regression with the ‘lm‘ function in R to assess the correlation between IQ and the number of ultra-rare class I variants in all autosomal genes and constrained genes (pLI ≥ 0.9). The analysis was adjusted for sex, first ten PCs from PCA, number of ultra-rare synonymous variants and total number of variants.

Overlap with common-variant risk loci

Common-variant gene-based associations were calculated using MAGMA⁷² and summary statistics from our previous GWAS meta-analysis of ADHD¹³. Association was tested using the SNP-wise mean model and linkage disequilibrium correction was based on estimates from samples with European ancestry from phase 3 of the 1000 Genomes Project⁷³. No window around genes was used. Gene-based results were subsequently used in a MAGMA competitive gene-set analysis to test for the enrichment of common-variant associations in two gene sets: (1) rare-variant risk genes with P < 0.001 (20 genes); and (2) rare-variant risk genes with P < 0.005 (62 genes).

Joint effect of common and rare variants

For the individuals in iPSYCH, common-variant data are also available. These data were included in our previous GWAS meta-analysis of ADHD, which contains detailed information on data generation¹³. In short, the samples were genotyped using Illumina’s PsychChip (iPSYCH1 samples) or Illumina’s Global Screening Array (iPSYCH2 samples). QC, imputation and association analysis were done using the bioinformatics pipeline Ricopili⁷⁴. After stringent QC, individuals and variants were included according to the following parameters: subject call rate > 0.95, autosomal heterozygosity deviation (|Fhet| < 0.2), variant call rate > 0.98, difference in variant missingness between cases and controls < 0.02 and SNP Hardy–Weinberg equilibrium (HWE) (P > 10⁻⁶ in controls and P > 10⁻¹⁰ in cases). Imputation was done separately for iPSYCH1 and iPSYCH2 samples using EAGLE v.2.3.5⁷⁵ and Minimac3⁷⁶, and the Haplotype Reference Consortium⁷⁷ panel v.1.0 was used as reference.

The PGS was constructed by splitting the relatedness pruned GWAS dataset (25,895 with ADHD; 37,148 controls) in 50 random subsets of roughly even size. For each of these, a GWAS was run on the complimentary 49 subsets using 10 PCs as covariates, and the results were meta-analysed with PGC and deCODE ADHD GWAS summary statistics described elsewhere¹³. The resulting summary statistics were then used to generate PGSs in the index set using the SBayesR algorithm implemented in LDAK⁷⁸. Finally, the scores from the 50 subsets, which together cover the full dataset, were assembled into one dataset.

We binned the PGSs into pentiles, and ten dummy variables were generated identifying the individuals in each pentile bin who have at least one class I variant or none, respectively. Logistic regression of ADHD status on the nine dummy variables was performed to estimate the impact on ADHD risk in each PGS pentile using individuals in the first pentile with no class I variants in constrained genes as reference. This was done separately for individuals with no class I variants in constrained genes and for individuals with at least one class I variant in constrained genes. The regression was adjusted for individual birth year, total number of variants, number of rare synonymous variants, percentage of exome target covered at a read depth of at least 20, mean read depth at sites within the exome target passing VQSR, sequencing wave, and the first ten PCs.

C-alpha test

To assess whether the rare variants identified in ADHD and ASD come from the same underlying gene distribution, we performed C-alpha tests on the iPSYCH exome data, including both ADHD and ASD. This comparison was designed to be simultaneous. To start with, we identified individuals diagnosed with either ADHD or ASD by the end of 2016, along with a control group in the quality-controlled iPSYCH data (n = 28,448). Subsequently, we defined and classified the rare variants into class I, class II and synonymous categories.

We performed a C-alpha test between ADHD without ASD (ADHD only) and ASD without ADHD (ASD only), regardless of ID comorbidity. In addition, we stratified the samples by the presence or absence of ID to perform C-alpha tests between ADHD only and ASD only. A separate set of C-alpha tests compared the single disorders with the control group.

The C-alpha test³⁶, implemented in the AssotesteR R package (http://cran.r-project.org/web/packages/AssotesteR/index.html), was used to evaluate the similarities between ADHD and ASD, and between controls and either ADHD or ASD. Each pairwise comparison for class I, class II and synonymous variants underwent 10,000 permutations, allowing us to verify the asymptotic P value against the permutation-based P value.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.