Patients with EOS were initially recruited through the 2011–2013 Interregional Hospital Clinical Research Program (NCT01512641). The main goal of this program was to estimate the prevalence of EOS in a population of children in child psychiatric care or medico-educational structures (Dor-Nedonsel et al., 2020). In addition, an enrolment in the GenAuDiss protocol (NCT02565524) study was offered to in- and out-patients, to Child and Adolescent Psychiatry (CAP) centers, and to their first-degree relatives (Fernandez et al., 2018). Participants were included either directly at the study sites or after referral by child and adolescent psychiatrists of the Provence-Alpes-Côte d'Azur (PACA) region of south of France. At inclusion (V1), clinical assessments (both psychiatric and neurocognitive) were performed in patients as well as first-degree relatives (parents and siblings). The following quantitative and qualitative measures were assessed during the study: 1) Clinical parameters from medical history (including pregnancy and birth) and biographic parameters, types and dates of significant life events including trauma and environmental exposures (e. g., drugs and substances); 2) Clinical parameters from physical examination (e. g., body weight, BMI and arterial pressure); 3) Semistructured interviews to assess the main diagnosis at the time of inclusion (K-SADS-PL, ADI-R, and MINI); 4) Clinical heteroassessments with specific neurocognitive rating scales (WISC-V, TMTA, TMTB, and Verbal fluency) and specific psychiatric rating scales (PANSS and SANS); self-report questionnaires (Cloninger TCI 226 to assess personality trouble and Baron Cohen AQ to assess autistic traits). In Table 1, all the main phenotypes are recapitulated.

Blood samples of trios (patients and their parents) were collected by the staff of the University Children’s Hospital of Nice (France). Samples were anonymized and DNA from blood was extracted using the GenElute Mammalian Genomic DNA Miniprep Kit (SIGMA Cat. No. G1N70). First, we performed current testing (including DNA testing for Fragile-X Syndrome, high-resolution karyotype and CGH-array). Then, in case of negative results from this step, we performed WES on trios. Isolated DNA samples were sent for WES analysis to IntegraGen SA (Evry, France).

Library preparation, exome capture, sequencing and data analysis were done by IntegraGen SA (Evry, France). Genomic DNA was captured using Agilent in-solution enrichment methodology (SureSelect SureSelect XT Clinical Reasearch Exome, Agilent) with their biotinylated oligonucleotides probes library (SureSelect XT Clinical Reasearch Exome—54 Mb, Agilent), followed by paired-end 75 bases massively parallel sequencing on an Illumina HiSeq4000 instrument.

Sequence capture, enrichment and elution were performed according to the manufacturer’s instruction (SureSelect, Agilent) without modification except for library preparation performed with the NEBNext^® Ultra kit (New England Biolabs^®). For library preparation 600 ng of each genomic DNA were fragmented by sonication and purified to yield fragments of 150–200 bp. Paired-end adaptor oligonucleotides from the NEB kit were ligated on repaired, tailed fragments then purified and enriched by 8 PCR cycles. 1,200 ng of these purified libraries were then hybridized to the SureSelect oligo probe capture library for 72 h. After hybridization, washing and elution, the eluted fraction was PCR-amplified with 9 cycles, purified and quantified by qPCR to obtain sufficient DNA template for downstream applications. Each eluted-enriched DNA sample was then sequenced on an Illumina HiSeq4000 as paired-end 75 base reads. Image analysis and base calling were performed using Illumina Real Time Analysis (2.7.3) with default parameters.

Base calling was performed using the Real-Time Analysis software sequence pipeline (2.7.3) with default parameters. Sequence reads were mapped to the human genome build (hg19/GRCh37) using Elandv2e (Illumina, CASAVA1.8.2) allowing multiseed and gapped alignments. Duplicated reads (e.g., paired-end reads in which the insert DNA molecule showed identical start and end locations in the human genome) were removed.

CASAVA1.8.2 was used to call single-nucleotide variants (SNVs) and short insertions/deletions (max. size = 300 nt), taking into account all reads per position. SNVs and indels with Q (SNPs) < 10 and Q (Indel) < 20 or regions with low mappability (QVCutoff <90) were filtered out. The frequency with which single base differences are expected between two unrelated haplotypes (Theta parameter) is 0.01. This frequency was set to 0.001 for indels.

Variant annotation takes into account data available in dbSNP (dbSNP144), the 1,000 Genomes Project (phase1_release_v3.20101123), the Exome Variant Server (ESP6500SI-V2-SSA137), and the Exome Aggregation Consortium (ExAC r3.0) and from an in-house database (201 exomes whole exomes for SNVs and 130 exomes whole exomes for indels). Functional consequences of variants on genes, transcripts and protein sequence, as well as regulatory regions, are predicted by Variant Effect Predictor (VEP release 83) (stop, splicing, missense, synonymous … ), as well as by location of the variants (e.g., upstream of a transcript, in coding sequence, in non-coding RNA, in regulatory regions). Regarding missense changes, two bioinformatic predictions for pathogenicity were used: SIFT (sift5.2.2) and PolyPhen (2.2.2). Other informations like quality score, homozygote/heterozygote status, count of variant allele reads, the presence of the variant in the COSMIC database (version71) were reported.

To investigate genomic copy number aberrations (CNA) (e.g., copy number gains and copy number losses), we used the Bioconductor DNACopy package (DNAcopy 1.32.0) by comparing the normal DNA exome data to a reference sample pool. It implements the Cystathionin Beta-Synthase algorithm to segment DNA copy number data. All changes were annotated with the catalog of the Database of Genomic Variants (DGV) to provide a comprehensive summary of structural variation in the human genome.

Gene prioritization was performed by combining the multiple bioinformatic tools we describe here. RVIS is based on the assumption that genes with common genetic functional variations are more tolerant and the variations occurring in these genes are common (higher RVIS score) and are unlikely to be disease-causing. On the other hand, genes with lower genetic variation (lower score) are considered to be intolerant and have high chances to be disease-causing. The RVIS percentile refers to the percentage of intolerant genes among which the gene of interest is found (Petrovski et al., 2015). GDI predicts if a given gene is prone to harbor disease-causing variants. The index evolved from the observation that more than a half of rare variants found in protein coding sequences are located in 2% of the genes (Itan et al., 2015). STRING is a protein-protein interaction database that selects data based on experiments (co-crystallization, co-purification and genetic interaction), databases (pathways and protein complexes), textmining (names of proteins that are often mentioned together), co-occurrence of gene patterns across genomes, and coexpression (Franceschini et al., 2013). VarElect is a tool that enables gene prioritization according to the disease and/or phenotype association. It also includes access to tools such as RVIS and GDI (Stelzer et al., 2016). DECIPHER is a database that combines multiple tools to facilitate evaluation of the impact of a given variant (Chatzimichali et al., 2015). The PANTHER classification system facilitates protein and gene categorization according to family and subfamily classes, molecular function, biological process, and pathways (Mi et al., 2020). MutationTaster enables to estimate the pathogenic potential of DNA sequence variants (Schwarz et al., 2014).

The two tails Mann-Whitney Test was performed for all the analyses. Data are presented in the form of violin plots. Statistical analyses were performed using the Prism Software 7 and 8 versions (GraphPad Software, Inc.).

The study protocol was approved by the Local Ethic Committee “Sud Méditerranée V” (number 14.002) and authorized by the French National Agency for Medicines and Health Products Safety (ANSM 2013-A01699-36). All patients signed a document authorizing the analyses described in this study. Collection of patients was authorized for a precise period (April 2014 - November 2021) and it was registered on ClinicalTrials.gov (NCT02565524).

As described in the Material and Methods section, we carried out WES on patients belonging to nine unrelated families whose genealogic trees are described in Figure 1.

Previous data suggest that there is a difference between the RVIS and GDI scores of inherited and de novo mutations. De novo mutations are thought to occur more frequently in genes that are more prone to harbor disease-causing variants and have less common variation (Ambalavanan et al., 2016). We compared the RVIS and GDI scores of inherited and de novo mutations, but statistical analysis with the two tails Mann-Whitney Test did not reveal statistical differences between these two groups (Figure 2). Even when we considered only missense de novo mutations (Supplementary Table S1), we did not obtain the same pattern as in previous studies (Figure 2). Actually, more missense variants occurred in genes that were predicted to have more common variation (Supplementary Table S1), which is opposite to what was found by other groups (Ambalavanan et al., 2016).

In order to find a link between trios, we created a Venn Diagram with the de novo variants identified in families 1–9 (Supplementary Figure S1). We focused on exon and splicing intron sequences, taking into consideration missense, nonsense plus splice exon, frameshift, in frame deletion/insertion, and splice acceptor variants. Due to splice variant identification, some changes were localized in the intronic part of the genes. We focused on de novo variants in genes that were common for at least three families from our dataset (Supplementary Table S2).

In families 1-5 and 7, we identified many de novo changes in Zinc Finger Protein 717 (ZNF717) that is a transcriptional regulator coding the KRuppel-Associated Box (KRAB) zinc-finger protein. According to the DECIPHER database, ZNF717 was found within the CNV region in five patients suffering from ID and in two patients with ASD. All reference SNPs (rs2918517, rs3009004, rs78135954, rs80214016, rs74357986, rs77110669, rs79138891, rs73843014, rs75737034, rs150497643, rs149076283 rs77378861, and rs776210532) do not have ClinVar reference. Most altered alleles have low frequency or data are not available. GDI predicts that changes within the ZNF717 gene have moderate chances to be disease-causing. In 2018, a de novo variant in ZNF717 was reported as a mutation potentially associated with growth retardation and ID (Diao et al., 2018).

Variants in Mucin 6 (MUC6) were identified in families 1–5. MUC6 is actually one of the most frequently mutated genes in the general population and are susceptible to false positive results in next generation sequencing http://massgenomics.org/2013/06/ngs-false-positives.html (Itan et al., 2015). However, evidence of late-onset Alzheimer disease (LOAD)-associated genetic polymorphism within an exon of Mucin 6 (MUC6) and immediately downstream from the Adaptor Related Protein Complex 2 Subunit Alpha 2 (AP2A2) have been reported (Katsumata et al., 2020).

Variants of the FSHD Region Gene 1 (FRG1) gene were found in Families 1–3. FRG1 (located in 4q35) codes a splicing factor and was implicated in FascioScapuloHumeral muscular Dystrophy (FSHD) as a modifier gene (DeSimone et al., 2017) (Jagannathan et al., 2019). 4p/4q CNVs of variable sizes have been associated with several different psychiatric findings and developmental disability. It was shown that the prevalence of independent deletions at 4p/4q involves PIGG, TRIML2, and FRG1. This suggests a possible implication of FRG1 in neurodevelopmental disorders (Novo-Filho et al., 2016). This gene is also susceptible of false positive results (http://massgenomics.org/2013/06/ngs-false-positives.html).

Another variant shared by families 1–3 is Poly(A) Binding Protein Cytoplasmic 1 (PABPC1). This gene was shown to be upregulated in male patients affected by SCZ (Qin et al., 2016). According to GDI prediction, PABPC1 has moderate chances to harbor disease-causing mutations. On the other hand, RVIS predicts that this gene has less common genetic variation.

Variants in Aggrecan (ACAN) were found in families 2, 3, and 4. The variant found in family 4 (rs12899191) turned out to have a high frequency in the general population [G = 0.27107 (3,538/13,052, GnomAD)] and is not reported in the ClinVar database. DECIPHER analysis shows that this gene is extremely intolerant to loss-of-function variation (Burd and Kerbeshian, 1987) (Gochman et al., 2011). ACAN was shown to be downregulated in the brain cortex region of SCZ individuals (Pietersen et al., 2014). Both RVIS and GDI predict this gene as highly mutated which means that there is little probability to be a disease-causing gene. A same kind of variant in RNA-Binding Motif Protein, X-Linked-Like-3 (RBMXL3), an RNA binding protein with elevated recurrent mutation rate in neuroblastoma (Lee et al., 2020), was found in families 1, 3, and 4. The DECIPHER database provides reports of two patients with abnormalities in the nervous system who carry variants in the RBMXL3 gene. MutationTaster shows that deletions in this gene might be disease-causing. GDI places this gene among ones that have moderate chances to harbor disease-causing mutations. The RVIS database does not contain information about this gene.  In order to group the genes harboring de novo variants according to their association with certain diseases or phenotypes, we carried out analysis with the VarElect tool. As a query, we used the following phenotypes: SCZ, ASD, ID, and neurodevelopmental disorders. We concentrated on the first 10 best scored genes (Table 2).

A missense mutation in the Potassium Voltage-Gated Channel Subfamily B Member 1 (KCNB1) gene was found in a patient of family 6. This mutation is predicted by the IntegraGen tool to have a moderate impact. SIFT and PolyPhen predict this change as deleterious (score 0) and probably damaging (score 0.999), respectively. MutationTaster indicates this change as disease-causing. KCNB1 is mainly linked to Early Infantile Epileptic Encephalopathy (de Kovel et al., 2017). The proband was, apart from EOS, diagnosed with idiopathic epilepsy, thus, this mutation identified in KCNB1 may be associated with the epileptic phenotype of the patient.

The second gene predicted by VarElect is Major Histocompatibility Complex, Class II, DR Beta 1 (HLA-DRB1). A missense variant (rs201929247) in this gene was identified in family 4. MutationTaster predicts this variant to be of polymorphic nature. Moreover, this variant turns out to be very frequent [A = 0.147,043 (14,841/100,930, ExAC)] and is not reported in ClinVar. According to GDI and RVIS predictions, variations in this gene are very common. For these reasons, its impact on the phenotype is very dubious.

A variant for the third predicted gene Bromodomain PHD Finger Transcription Factor (BPTF) was identified in family 8. BPTF is associated with a neurodevelopmental disorder characterized by ID, dysmorphic facies and distal limb anomalies (Midro et al., 2019). The IntegraGen tool predicts the impact of this change (in-frame insertion) as moderate. GDI analysis suggests that BPTF is among moderately-damaging genes. RVIS predicts this gene as a gene of less common variation.

A missense variant identified in HECT And RLD Domain Containing E3 Ubiquitin Protein ligase 2 (HERC2) was found in family 1. Mutations in HERC2 were associated with developmental delay with Angelman-like features (including severe ID) (Harlalka et al., 2013). This variant (rs146883683) is known to have high frequency [G = 0.11394 (9,547/83,792, ExAC)] and is not reported in the ClinVar database. Moreover, MutationTaster predicts this variant as polymorphism.

A missense variant in Cystathionin Beta-Synthase (CBS) was identified in family 3. According to ClinVar, the variant rs5742905 is associated with homocystinuria. Interestingly, this disorder is frequently accompanied by ID and CBS is one of the genes associated with Down Syndrome (Pogribna et al., 2001), (Marechal et al., 2019). MutationTaster predicts this change as disease-causing. GDI analysis suggests that CBS is among moderately damaging genes and RVIS predicts this gene as a gene of less common variation.

An insertion in Neurogenic Locus Notch Homolog Protein 4 (NOTCH4) was identified in family 4. MutationTaster predicts this variant as a polymorphism. Interestingly, polymorphic variants in this gene were previously associated with SCZ and schizoaffective disorder (Ujike et al., 2001), (Zhang et al., 2015). Both GDI and RVIS predictions indicate that variants in this gene are not prone to be disease-causing.

A missense variant in Folate Hydrolase 1 (FOLH1) was found in family 3. This variant (rs75111588) is frequent in the general population [A = 0.187,153 (20,786/111,064, ExAC)]. Hence, it is unlikely that this variant is responsible for the phenotype of the patient.

A missense variant in Filaggrin (FLG) was identified in family 5. It turned out to be frequent in the general population [G = 0.08408 (1,087/12,928, GO-ESP; G = 0.117 (70/600, NorthernSweden)], which suggests that this variant does not impact the phenotype. FLG mutations were associated with athopic dermatitis with increased ischemic stroke risk in the general population (Varbo et al., 2017).

An in-frame deletion in Achaete-Scute Family BHLH Transcription Factor 1 (ASCL1) was identified in family 6. The IntegraGen tool predicts this change to have a moderate impact. GDI places this gene among the genes with moderate chances to harbor disease-causing variants. ASCL1 promotes neural differentiation (Castro et al., 2011) and its level is more elevated during neural differentiation in murin Fmr1-null embryonal stem cells compared with WT cell line (Khalfallah et al., 2017). According to the genomAD browser this variant is quite frequent (delGCA = 0.01150; 340/29,556). MutationTaster predicts this change as potential polymorphism.

In family 6, we also identified a frameshift variant in Zinc Finger Protein 148 (ZNF148). IntegraGen predictions indicate high impact of this variant. Mutations in ZNF148 were associated with Global Developmental delay, Absent or hypoplastic Corpus Callosum and dysmorphic Facies (GDACCF). Apart from this disorder, the phosphorylation level of ZNF148 was shown to be decreased in SCZ patients, but its overall level of proteins was unchanged (Jaros et al., 2012).

In general, the VarElect analysis has provided a further insight into genes that are already implicated in the pathogenesis of neurological disorders. We could find interesting variants in genes such as ZNF148 and KCNB1 that could impact the phenotype of the proband of family 6. We were able to exclude some variants of polymorphic nature that turned out to be common in the general population and that were not associated with any neurological disorder. We came across polymorphisms that could increase the risk of SCZ (NOTCH4). The VarElect analysis enables to extract variants that sometimes can be overlooked in general analyses (e.g., some polymorphisms that are frequently excluded from candidate genes).

We created a Venn Diagram for inherited variants to find genes that are common among at least three families (Supplementary Figure S2). We found shared inherited variants in the following genes: ZNF717, FAM182B, LRP5L, HRNR, OR4C5, APOL3, BCLAF1, FLG, RP1L1, MUC6, HYDIN, PLEC, TTN, NRAP, ADAMTS8, SSPO, FAM104B, DAAM2, ZFHX3, and AHNAK2. According to GDI and RVIS predictions, most of these genes are less prone to harbor disease-causing mutations than the genes cited above (de novo variants) and belong to genes that have more variability.

We focused on genes that did not have RVIS references. In families 1 to 7, we found 17 missense variants, two frame-shifts and one insertion in ZNF717. This gene was previously found during de novo mutation analysis. Considering the high amount of variants found and the fact that the majority of them have references in SNPdb, variants in this gene may be very common and not prone to be disease-causing.

In families 3 and 4, we found variants in Family With Sequence Similarity 182 Member B (FAM182B). In family 4, we found a variant in the 5′-UTR and a missense variant. In family 3, we found a variant in the 5′-UTR and one synonymous variant. This gene is predicted by GDI to have medium chances to harbor disease-causing mutations.

In families 1, 2, and 9, we identified variants in Hornerin (HRNR). All of them have references in SNPdb (rs79513582, rs76694305, rs61814943, and rs61814936) and have high frequencies in the general population.

We found some variants in Olfactory Receptor Family 4 Subfamily C Member 5 (OR4C5): 1) the same deletion (rs66829866) in families 1 and 3 was reported in SNPdb and is known to have a high frequency; 2) an insertion in family 1 (rs114053360) that has low frequency and no report in ClinVar; 3) an insertion in family 4. OR4C5 is a gene predicted by GDI to be highly damaging.

We identified three variants in Plectin (PLEC). A missense variant in family 2 is predicted by PolyPhen2 to be benign (score 2). In family 4, we found a synonymous variant (rs113133985) described in SNPdb with low frequency. A missense variant in family 5 also turned out to have low frequency and no records in the ClinVar database. GDI was predicted to be high, whereas RVIS was very low.

Overall, all inherited variants that we found seem to have benign consequences.

We carried out interaction analysis with the STRING database (Supplementary Table S3). Due to the fact, that some mothers from our cohort have mild symptoms associated with neurodevelopmental disorders, we decided to focus on the genes inherited from the mothers of the probands. We took into consideration all inherited genes regardless of their potential impact of the carrying variants. The analysis of the mutated genes that were inherited from the mother revealed that “the network among the genes subjected to the analysis had significantly more interactions than expected” (p = 7.08 × 10⁻⁶). This means that these proteins have more interactions among themselves than what would be expected for a random set of proteins of similar size drawn from the genome. Such an enrichment indicates that the proteins are at least partially biologically connected. Molecular functional enrichment in this network includes transferase activity [False Discovery Rate (FDR) = 0.00019], catalytic activity (acting on protein, FDR = 0.0011), adenyl ribonucleotide binding (FDR = 0.0011), ATP-binding (FDR = 0.0011), and ubiquitin protein ligase activity (FDR = 0.0014). For instance, among the genes in the enrich pathways we came across POLA1, PARP1, and CDK5 (full gene list in Supplementary Table S1). Interestingly, POLA1 was recently associated with X-Linked ID associated with severe growth retardation, microcephaly, and hypogonadism (Van Esch et al., 2019). Elevated levels of PARP1 have been associated with Down Syndrome which is characterized by ID (Salemi et al., 2014) and it was involved in the context of dyslexia and neuronal migration (Tapia-Páez et al., 2008). Dysregulations of the CDK5 activator are associated with SCZ (Engmann et al., 2011). Considering PFAM protein domains, enrichment was observed in protein kinase (FDR = 0.0081), protein tyrosine kinase (FDR = 0.0089), microtubule-binding (FDR = 0.0439), and kinesin motor domains (FDR = 0.0439).

We performed pathway analysis with PANTHER for altered genes that were inherited from the mother. We looked for genes that belong to the same pathways. We assigned each gene to the corresponding pathway according to PANTHER classification. Subsequently, a Venn diagram was created in order to see if the pathways co-occur in our cohort of families (Supplementary Figure S3). Functional classification analysis showed that all nine families share variants in genes belonging to the WNT and Cadherin pathways (Supplementary Figure S3). In addition, families 1 to 8 share additionally the cholecystokinin receptor (CCKR) signaling pathway. Cholecystokinin activates the CCK1 receptor that belongs to G-protein-coupled receptors. Variants associated with these pathways by PANTHER are comprised, among others, in the following genes: CDH11, CDH3, CDH18, PCDH15, PCDHA4, FAT4, FAT2, NFATC2, WNT6, WNT11, and CTNNA3.