We used the DDD data acquired from 24 clinical genetics centers within the United Kingdom (UK) National Health Service and the Republic of Ireland. Together, they recruited 13,462 individuals based on undiagnosed developmental disorders and their parents (9,860 trios) from April 2011 through April 2015 (16, 17). The inclusion criteria for the DDD study to recruit children is based on the evidence of neurodevelopmental disorder, congenital anomalies, abnormal growth parameters (height, weight, occipitofrontal circumference), dysmorphic features, unusual behavioral phenotype or genetic disorder of significant impact for which the molecular basis is currently unknown (18). After obtaining ethics approval at our centers, and with permission from the DDD consortium, we analyzed the dataset for phenotypes of the individuals. We downloaded and worked on the GrCh37 genome-aligned files from the last datafreeze of 2017-12-15.

The phenotypes listed in DDD follow the Human Phenotype Ontology (HPO) organization. The hierarchical structure of HPO can be accessed in two formats: OBO flat file format and Web Ontology Language (OWL). The OWL version is more enriched as it provides logical definitions of the HPO classes, which facilitates cross-species mapping with semantic reasoning tools (19). OWL format was used as it consists of existential and universal restrictions on HPO phenotype classes. Protégé, which is an open source ontology editor and framework (20), was used to explore the defined restrictions on different phenotypes. These restrictions describe anonymous ancestors of the phenotypes and were used to classify the phenotypes according to the body section. HPO presents the advantage of being a machine readable and cross-species usable ontology but has limitations when being used to characterize fine clinical phenotypes.

We identified individuals with phenotypes that included global developmental delay, which were then divided by sex. We excluded people whose sex or age was not specified, so a number of 6,988 individuals with GDD are retained for this study.

The exome sequence data of GDD phenotyped individuals were analyzed in two stages. In the first stage, the existing GRCh37/hg19 exome sequences were realigned to the GRCh38 genome reference sequence. Then the short variants (SNVs and indels) were called using GATK best practices workflow (21), which involves steps such as realigning reads to GRCh38 reference genome, variant calling using HaplotypeCaller and joint genotyping, and finally variant quality recalibration and refinement steps leading to a high quality variant callset. In the second stage, these variants were annotated for gene information (Ensembl), frequencies (from gnomAD, ExAC, and internal cohort GDD), and pathogenicity (from CADD, Clinvar, and Clingen). The annotated set of variants in the callset were filtered for gene information, rare variants having minor allele frequency (MAF) <=0.01, impact on the transcript, and pathogenicity (pathogenic/likely_pathogenic in CLINVAR/CLINSEG). The details of the annotation and filtering criteria can be found in the Supplementary data Section S1.1. For some analyses we considered genes with a higher probability of loss-of-function intolerance score (pLI >= 0.9).

In order to assess genomic data in a targeted manner, we developed disorder specific candidate gene lists. We searched PubMed using the keywords: intellectual disability (ID)” and “global developmental delay (GDD)”, reviewed the papers and compiled manually a list of genes from original research and review papers (22–25). We integrated both GDD and ID together as most papers reviewed did not make a clear distinction between the two. We also used genes included in already developed in trusted databases related to neurodevelopmental disorders (NDD) for diseases, phenotypes, and even genome wide association studies (GWAS) using the same keywords [SysID (26, 27), DisGenet (28, 29), HPO (19, 30), OMIM (31, 32), Orphanet (33), Phenolyzer (34, 35), Ingenuity Pathway Analysis (Qiagen), Open Targets (36, 37), AutDB (38)]. We also added the Intellectual Disability NGS Radboudumc and Fulgent gene panels to have the most complete overview. Each gene list was obtained separately, and then we retained only those genes that appeared at least 3 times in the collected data, resulting in a list of 2,539 candidate genes for ID/GDD (Supplementary Table S1a). We also similarly created a list of autism spectrum disorder (ASD) candidate genes. We searched PubMed using the keywords: “autism”, “autism spectrum disorder ASD”, and “autistic” and compiled a list of genes from original research and review papers (23, 39). To the database list used for ID + GDD gene list we added the database created by SFARI (40), resulting in a list of 730 genes (Supplementary Table S1b). We mainly based ourselves on curated databases from autism experts (SFARI) considering the issues with annotation which has been seen when using original reports (41). Several articles have already been screened by experts to be focused on ASD as some papers present genes related to autism but also associated with other traits. We are also using a stringent approach by using a threshold of 3 references per gene.

Expression data was retrieved from GTEx v8 brain tissues to compare expression in males and females in brain regions (11 distinct brain regions and 2 cell lines) (42), using the medium Transcripts Per Million (TPM). We also used the BrainSpan data from both sexes, containing 19 brain regions (43). The median Reads Per Kilobase per Million (RPKM) was reported for each sex by tissue before and after birth.

We used logistic regression to highlight the differences between males and females with respect to phenotypes. Briefly, the logistic regression model estimates odds-ratios for binary outcomes. Only phenotypes with N > 5 in each category were retained in the analysis. Since we expected possible overlaps between phenotypes, due either to common causal mechanisms or phenotype misclassifications, we applied the Benjamini Yekutieli correction to adjust p-values with the false discovery rate (FDR) for the multiple comparisons. Adjusted p-values lower or equal to 0.05 were considered statistically significant. We next used a linear mixed model to take into account the age in the model for the adjustment of the difference of the phenotypes between the sexes.

Consistently, to highlight the differences between males and females regarding the mutations in candidate genes, we still used a logistic regression model. However, due to strong mutated gene imbalance between sexes, we observed a bias toward conservatism. To address this limitation, nominal p-values lower or equal to 0.05 were considered statistically significant.

To assess whether the gene expression is impacted either by sex or developmental stage, we took advantage of linear mixed model of the logarithm of gene expression conditioning on sex (reference = female) and sex specificity (reference = male specific) for GTEx data, and for sex (reference = female), sex specificity (reference = male specific) and developmental stage (reference = after-birth), when analyzing BrainSpan data. One advantage of linear mixed models over others approaches is to model the dependence structure of expression for the same gene, while removing unmeasured confounding factors. Thus, we considered a random intercept for each gene. It is worth noting that for these analyses we adjusted for multiplicity using Bonferroni correction. Adjusted p-values lower or equal to 0.05 were considered statistically significant.

Finally, cluster analysis was performed using hierarchical cluster analysis methods with Ward distance. We chose this metric since it minimizes the intra-cluster variance while maximizing the inter-cluster variance, giving clusters with genes of closed expression. Threshold was chosen in both males and females for each brain region, using a graphical criteria.

For each group (males or females with GDD or other phenotypes of interest) the genes harboring mutations were then analyzed by protein-protein interaction prediction, clustering, and pathway enrichment using STRING v11.5 with full STRING network, all active interaction sources and medium confidence (0.4) (44). The Markov Clustering Algorithm (MCL) was used to identify gene subnetworks (45). The inflation parameter of MCL was set to 1.5 (46). Functional enrichment analysis of each module was also performed in the STRING v11.5 database using the GO terms, REACTOME pathways, and KEGG pathways. The false discovery rate obtained from the functional enrichment analysis describes the degree of significance of the enrichment. The p-values were corrected for multiple testing in each category using the Benjamini-Hochberg procedure. Visualizations of protein-protein interactions and clusters were also obtained using STRING v11.5.

In order to capture biological pathways occurring in male specific and female specific gene clusters, we applied cluster-based gene-ontology enrichment analyses using the function compareCluster from clusterprofiler R package (47).

We observed significant differences between sex in several phenotypes in individuals with GDD (Figure 1, Supplementary Table S3). Indeed, we observed that males with GDD had a significantly higher rate of autistic behavior (HP:0000729, 7.57% males, 4.39% females, P = 9.82E-05) and autism (HP:0000717, 5.62% males, 3.09% females, P = 4.22E-04). We also found a significant increase in stereotypical movements (HP:0000733, 3.45% females, 1.54% males, P = 3.31E-04), short philtrum (HP:0000322, 2.87% females, 1.51% males, P = 4.70E-02), and congenital hip dislocation (HP:0001374, 0.94% females, 0.18% males, P = 7.15E-03) in females. There was a significant increase in frontal upsweep of hair (HP:0002236, 1.56% males, 0.40% females, P = 1.13E-03) and inguinal hernia (HP:0000023, 2.45% males, 0.25% females, P = 1.37E-11) in males. We also used logistic regression to assess whether there was a difference when adjusting for sex and age. We observed no significant differences after adjusting for multiple comparisons. Also, the p values of the interaction of the age through the sex were not significant. This means that there is insufficient evidence to conclude that the interaction of age through sex influences the phenotypes seen above (Supplementary Table S5).

Interestingly, we found several phenotypes where males and females presented with quantitatively or qualitatively “opposite” phenotypes (Figure 1). For example, down slanted palpebral fissure was increased in males (HP:0000494, 5.44% males, 3.23% females, P = 7.08E-03) whereas upslanted was increased in females (HP:0000582, 5.88% females, 3.25% males, P = 2.54E-04). Also, microcephaly was significantly more frequent in females (HP:0000252, 18.22% females, 12.88% males, P = 4.58E-06), while macrocephaly was significantly more common in males (HP:0000256, 4.29% males, 2.58% females, P = 4.70E-02). Developmental regression was significantly increased in females, probably due to the inclusion of individuals with Rett syndrome which is found most commonly in females (HP:0002376, 2.32% females, 1.07% males, P = 2.38E-022.38E-02).

The rationale for DDD was to understand the genetic basis of undiagnosed individuals with NDD with more advanced genetic testing. Some individuals may have been included based on lack of clinical diagnosis (phenotypic presentation not fitting classical syndrome potentially) and then found to have mutations in genes syndromic. We did not exclude any condition, but the way in which the individuals in this database were recruited led us to the conclusion that there were no FXS in the cohort, as these were individuals with severe developmental delays who had not been diagnosed with any known disease. This could partly explain why we didn't observe a significant difference in the prevalence of GDD between the sexes, as FMR1 is the most common single gene causing developmental differences, which is biased in favor of males. Indeed, we didn't identify any individuals carrying a pathogenic variant in the FMR1 gene in the cohort. On the other hand, we found variants in the MECP2 gene, which corresponds to Rett syndrome, although we believe that individuals with this syndrome were excluded from the initial cohort. We analyzed further differences in genetic mutations between sexes to gain a better understanding of differences in phenotypes. We identified likely pathogenic and pathogenic variants in the individuals with GDD using ClinVar (49). We found that 66.3% of GDD genes, affected by pathogenic or likely pathogenic variants, were found in both males and females (Figure 2A), but also identified genes (with pathogenic or likely pathogenic variants) specific to males (20.4%) and females (13.2%). We assessed if the pattern of male and female specific affected genes corresponded to X-linkage and found that 8.3% of male specific genes were X-linked compared to 7.5% in females (Figure 2B).

We wanted to assess if differences in mutation load were present in GDD as previously reported in ASD (50). We therefore quantified the number of pathogenic mutations per individual in each sex (Figure 3A). We observed a similar number of pathogenic variants in GDD/ID genes per individual in males and females (median = 3, P = 0.84). We also considered GDD/ID genes present in both sexes but biased toward one sex (Supplementary Table S4). None of the genes that are significantly different between the two sexes are X-linked. Finally, we analyzed the genes having a higher probability (>0.9) of loss-of-function (LoF) intolerance (pLI) score. High pLI genes have been shown to have a higher impact on intellectual quotient in copy number variants in individuals with developmental differences (51, 52). Interestingly, we found that for those with high pLI genes, females had a significantly higher number of pathogenic variants per individual (0.59 vs. 0.52 for male, p = 1.7E-4) (Figure 3B).

Next, we hypothesized that differences in phenotypes could be related to differences in gene networks and molecular pathways. GDD genes carrying mutations present in both males and females subdivided in 27 sub-networks and 12 with more than 10 genes (Figure 4A). The first 3 sub-clusters, with the largest number of genes are shown in Figures 4B–D, which were involved in protein binding, chromatin binding and transcription factor; RNA polymerase, DNA binding; and catalytic activity in mitochondria (Supplementary Table S6). On the other hand, GDD genes found mutated specifically in males (Figure 5) also subdivided in 20 sub-networks including 4 with more than 10 genes. In the same way the first 3 sub-clusters were involved in DNA binding; gated channel activity; catalytic activity in mitochondria (Supplementary Table S7). Genes found mutated specifically in females (Figure 6), subdivided in 18 groups and the first 5 with more than 10 genes were associated with protein binding; gated channel activity; Endoplasmic reticulum activity; and mitochondrial enzymatic activity. So there were no clear cut differences in type of molecular function or cellular component between genes specific to males vs. females. (Supplementary Table S8).

We hypothesized that phenotypic differences in individuals with GDD may be related to differences in sex-biased gene expression in the brain, as was postulated before for ASD (53, 54). We leveraged sex-specificed gene expression data from two datasets: GTEx (42) and BrainSpan (43).

Firstly, we wanted to assess if regional differences in gene expression could explain phenotypic differences in GDD as seen recently in ASD (55). We assessed the gene expression profile by sex in each brain region, while adjusting for the GDD genes found mutated in males vs. females specific (and enriched) in the individuals with GDD (Supplementary Figure S1). We took advantage of the GTEx data in 11 brain regions (Amygdala, Caudate, Cerebellar Hemisphere, Cerebellum, Cortex, Hippocampus, Hypothalamus, Nucleus accumbens, Putamen, Spinal cord, and Substantia nigra) and 2 cell lines (Frontal Cortex (BA9) and Anterior cingulate cortex (BA24)). Using a linear mixed model of the logarithm of gene expression (See Methods), we found significant associations for the sex in all the considered brain regions. Moreover, when assessing the effect direction we observed that only in the amygdala, the anterior cingulate cortex, and the nucleus accumbens, being a male significantly decreases the gene expression, while for the remaining brain regions the association is positive. These results were globally confirmed when considering only sex specific genes (Results not shown). The results for all genes were reported in Supplementary Table S8. Then, we constructed gene expression clusters for sex specific genes and found different biological processes involved between the two sexes (See Methods). For example, focusing on the amygdala, the cerebellar hemisphere and the hippocampus, brain regions known to be strongly associated in emotions, higher cognitive functions, learning and memory, respectively, we found a weak overlap in GO pathways enrichments between males and females (Figure 7; Supplementary Figures S2-S7). Collectively, we found that sex is associated with gene expression across a large set of brain regions pointing to different biological processes between males and females. However, the crossed effect of sex and developmental stage on gene expression remains to be understood.

Another way genes can influence phenotypic presentation is via timing of expression in brain development as seen before in ASD (56). In order to address the role of the crossed effect of sex and the developmental stage on gene expression, we leveraged information from the BrainSpan database (43) for 19 brain regions [amygdaloid complex, anterior (rostral) cingulate (medial prefrontal) cortex, cerebellar cortex, dorsolateral prefrontal cortex, hippocampus (hippocampal formation), inferolateral temporal cortex (area TEv, area 20), mediodorsal nucleus of thalamus, orbital frontal cortex, posterior (caudal) superior temporal cortex (area 22c), posteroventral (inferior) parietal cortex, primary auditory cortex (core), primary motor cortex (area M1, area 4), primary somatosensory cortex (area S1, areas 3,1,2), primary visual cortex (striate cortex, area V1/17), striatum, ventrolateral prefrontal cortex, cerebellum, dorsal thalamus, and primary motor-sensory cortex (samples)]. We integrated the developmental stage in our model, while adding an interaction term between the developmental stage and sex to highlight synergies between these two variables on gene expression patterns. It is worth noting that we removed the dorsal thalamus, cerebellum, and primary motor-sensory cortex (samples), since we did not have complete data for our analysis. Firstly, we found that developmental stage is positively associated with gene expression in 4 brain regions: the amygdaloid complex, the cerebellar cortex, the striatum, and the hippocampus. These 4 regions have been shown to be strongly involved in a large range of developmental disorders (57–61). Moreover, sex-developmental stage interaction is significantly positively associated in the dorsolateral prefrontal cortex and the posterior (caudal) superior temporal cortex (area 22c) brain regions. Results were recapitulated within Supplementary Table S10. Consistent with results found in the GTEx data, these findings suggest that developmental timing may vary across sexes, resulting in different gene expression patterns among males and females, leading to phenotypic variations between sexes.

Next, we assessed a subgroup of individuals with GDD also presenting with autism, where sex differences have been observed (62). Autism is significantly overrepresented in males with GDD (5.62% in males vs. 3.09% in females P = 4.22E-04) (Figure 8A), but not in the 3–4:1 (62) found in the previous reports investigating individuals with ASD. First, we assessed if autism was more prevalent in males with GDD due to increased ASD gene mutations in males. We found that, actually, females had a higher mutation rate in ASD genes than males (0,74 in females vs. 0,68 in males; p = 0.0043, T-test) (Figures 8B,C). We then assessed if different types of ASD genes were mutated in males vs. females with GDD and autism. We found that 28 ASD genes were mutated in both sexes while 47 and 17 were unique to males and females, respectively (Figure 8D). We then examined their gene network (Figure 8D). The ASD genes shared by males and females were subdivided into 3 subnetworks with only one with more than 10 genes. Subcluster1 is related to chromatin binding, and it may be noted that one of the biological processes is memory (Supplementary Table S11). The ASD genes that are specifically mutated in males are divided into two subnetworks, one of which is larger and enriched for genes involved in voltage-gated ion channel activity and regulation of cell signaling, transport and activity (Supplementary Table S11). We can also note processes related to the synapse, neuron and dendrite. On the other hand, female specific genes did not show as much interaction between genes with only 2 small sub-networks with low functional enrichment which is related to methylation (Supplementary Table S11).

Next, we wondered whether sex could influence phenotypic presentation for the same gene. We tested if the phenotypic presentation of individuals with GDD was found to have pathogenic variants in ASD candidate genes (Supplementary Table S2) differed by sex. Interestingly, we found that males presented with autism, whereas females presented with an increase in stereotypies (Figure 8E). Similar results to those described above were obtained for genes with a high pLI (Supplementary Figure S8).