A total of two autistic children (IV2 and IV4) in a family pedigree were recruited between January 2019 and December 2020 for the present study (Fig. 1). They were diagnosed at the Second Affiliated Hospital of Wenzhou Medical University (Wenzhou, China) and identified as ASD according to Wenzhou Medical Identification Committee of Sick and Disabled Children. The diagnosis of ASD was confirmed by the Autism Diagnostic Observation Schedule (ADOS) (19) and expert clinical judgment. IV2 and IV4 children were found to have ASD characteristics at age 2 and 3 and diagnosed as autistic children at age 3 and 4, respectively. The scores of ADOS (model 1) were 10 and 8, respectively. IV2 was associated with mild cerebral palsy and mild intellectual disability, whereas IV4 was associated with epilepsy, brain development delay and mild intellectual disability. The parents of IV2 and IV4 were without clinical autism symptoms and the marriage type was not consanguineous. Psychological training, such as training in language communication and social interaction, was performed at the Mental Health Center, Wenzhou Medical University (Wenzhou, China), which is one of the author-affiliated institutes associated with this study. The present study was approved by the Ethics Committee of Wenzhou Hospital of Integrated Traditional Chinese and Western Medicine (approval number: 2018-57; Wenzhou, China); the Ethics Committee of this hospital was shared with the Wenzhou Maternal and Child Health Guidance Center (Wenzhou, China). Written informed consent forms was obtained from the parents of the children, for both their participation and the participation of the children. The clinical data and 2–3 ml/each peripheral blood samples of two autistic children and their parents were provided by Wenzhou Maternal and Child Health Guidance Center. The peripheral blood samples of IV2 and IV4 autistic children were collected at the age of 7 and 5 years, respectively.

WES of the autism samples (2 ml) was performed by iGeneTech Biotech (Beijing) Co. Ltd. Genomic DNA was extracted using a Blood Genomic DNA Isolation kit (magnetic bead method; Kangwei Biologic Co. Ltd.) in accordance with the manufacturer's protocols. The quality and quantity of DNA extracts were evaluated. DNA quality was examined using 1% agarose gel electrophoresis. The DNA concentrations were determined using a Qubit dsDNA BR assay kit (Thermo Fisher Scientific, Inc.). The whole exome liquid-phase capture technique was used to efficiently and specifically enrich all the human exonic regions. Briefly, fragmentation was carried out using a hydrodynamic shearing system (Covaris, Inc.) to generate 150–200 bp fragments. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities using the Fast Library Prep Kit [iGeneTech Biotech (Beijing) Co. Ltd.]. After adenylation of 3′ ends of DNA fragments, adapter oligonucleotides were ligated using the Fast Library Prep Kit. The Fast Library Prep Kit was used in accordance with the manufacturer's protocols. DNA fragments with ligated adapter molecules on both ends were selectively enriched by PCR. Subsequently, libraries were hybridized in liquid phase with a biotin-labeled probe, and magnetic beads with a streptomycin monolayer were employed to capture the exons of genes using the TargetSeq One Kit [iGeneTech Biotech (Beijing) Co. Ltd.], in accordance with the manufacturer's protocols. Captured libraries were enriched by PCR to prepare for sequencing. DNA library construction was performed using a Fast Library Preparation kit [iGeneTech Biotech (Beijing) Co. Ltd.]. Library quantification and quality assessment were then performed. Qubit 3.0 (Thermo Fisher Scientific, Inc.) was used for library quantification; library concentration >25 ng/µl was considered a qualified library. DNA quality was examined using Bioptic Qsep 100 DNA fragment analyzer (BiOptic, Inc.). The length of library fragments was generally between 270 and 350 bp. The Illumina NovaSeq6000 platform (Illumina, Inc.) was used for PE150 high-throughput double-terminal sequencing. The loading concentration of the final library was 450–600 pmol/l; library quantification was performed by qPCR. The raw WES data were obtained and mapped to the human reference genome (UCSC GRCh37/hg19, 2009). The type of analysis was single sample whole-exome analysis, which included quality evaluation of sequencing data as well as the detection and screening of variations. Based on the BAM results of alignment with the genome reference sequences, Samtools (http://www.htslib.org/) and GATK (https://gatk.broadinstitute.org/hc/en-us) were used to find single nucleotide polymorphisms (SNPs) and small insertions and deletions (InDel) of <50 bp. ANNOVAR (http://www.openbioinformatics.org/annovar/) software was used to annotate the SNP and InDel loci to determine the corresponding gene information, functional data and harmfulness. The annotations on SNP and InDel were carried out by iGeneTech Biotech (Beijing) Co. Ltd.

A large number of genetic variants were detected. To identify the variants with significant functions, the following conditions were set: i) The variants could be divided into high level, possible high level, intermediate level and low level according to the priority of harmfulness. In the present study, the screening for the priority of variants was limited to high level only. The high priority of the variant indicated that it was located in the exonic region of the gene, had a minor allele frequency (MAF) of <0.01 in East Asian populations [i.e. the 1000 Genomes Project database (1000G), the Exome Aggregation Consortium (ExAC) and the Genome Aggregation Database (gnomAD)] (20) and could lead to an altered amino acid sequence. In addition, at least 2 of the 10 prediction software packages [including functional prediction software SIFT, Polyphen2_HDIV, Polyphen2_HVAR, LRT, MutationTaster, MutationAssessor, FATHMM, fathmm-MKL, total effect prediction software CADD, conservative prediction software GERP; all from dbNSFP version 3.5a (annotated using Annovar; updated on September 21, 2018)] predicted it as harmful. ii) The quality evaluation of the variants could be divided into high, medium and low quality. In the present study, the quality of the variant was limited to high or medium quality only. The high quality of the variant is denoted by its depth of ≥50 layers, genotype quality value of >30 and frequency deviation of <0.1 (heterozygous 0.4–0.6; homozygous 0.9–1.0) (20). The medium quality of the variant is denoted by its depth of ~20 layers, genotype quality value of >20 and frequency deviation of <0.2 (heterozygous 0.3–0.7; homozygous 0.8–1.0) (20). iii) The mutant genes were required to be associated with autism in the human gene mutation database (HGMD) (http://www.hgmd.cf.ac.uk/ac/index.php).

All the preliminary screening variants of IV2 and IV4 were verified by Sanger sequencing. In addition, the parents of IV2 and IV4 were sequenced for these variants to determine the variant source. Conventional PCR amplification and Sanger sequencing were performed on all the exons where the variant loci of the corresponding gene were located; PCR primer synthesis and Sanger sequencing were conducted by Qingke Biology Co., Ltd. The primer sequences are provided in Table SI.

The AutDB autism database (http://www.mindspec.org/autdb.html; also known as SFARI Gene, http://gene.sfari.org) is updated in real time (21). By September 2020, it contained 1,237 genes related to autism, with an evidence score of 0–5. In the present study, the genes with evidence score of ≥3 in the AutDB database were selected as candidate autism-related genes. AutismKB Database is the world's largest database for the genetic evidence of ASD (3,4), with 1,379 genes and 11,669 single nucleotide variants or small InDels integrated in version 2.0 (2018). According to AutismKB, in addition to 99 syndrome genes, 1,280 non-syndrome genes with an evidence score of >4 were designated as candidate autism-related genes, of which 30 syndrome genes and 198 non-syndrome genes with an evidence score of >16 were designated as substantially credible autism-related genes. The criteria set by Yang et al (3) used for calculating the ‘evidence score’ of AutismKB v2.0 database were as follows: Next-generation sequencing (NGS) de novo mutation studies or low-scale gene studies, weight 10; NGS other studies, weight 8; genome-wide association studies, weight 4; Copy-number variation/structural variation studies, linkage analyses or low-scale genetic association studies, weight 2; and expression profiling or NGS mosaic mutation studies, weight 1.

Among all the detected variants in autistic child IV2, there were 21 variants (in 21 genes) with high harmful effects and high or medium quality, and their mutant genes were associated with autism in HGMD; further details are provided in Table I. Notably, the variant types of all genes were missense (nonsynonymous), except that GPR152 and OAS3 genes contained deletion (frameshift, InDel below 50 bp) and nonsense (stopgain) variants. All the variants were located in exons. The genotypes were all heterozygous and the variant frequency of each gene in this sample ranged from 0.32–0.59; that is, the minimum variant frequency (that of gene NRP2) was 0.32 and the maximum variant frequency (that of gene OAS3) was 0.59 in this sample. All of the mutant genes were located on autosomal chromosomes, except that ARHGEF6 was on the X chromosome.

Among all the detected variants in autistic child IV4, there were 22 variants (in 21 genes) with high harmful effects and high or medium quality, and their mutant genes were associated with autism in HGMD; further details are given in Table II. All of the variants were missense type and were located in exonic regions. The genotypes were all heterozygous, and the variant frequency of each gene in this sample ranged from 0.32–0.57; that is, the minimum variant frequency (that of gene MTMR12) was 0.32 and the maximum variant frequency (that of gene SMARCC2) was 0.57 in this sample. The mutant genes were all located on autosomal chromosomes. Based on the comparative assessment of WES results in Tables I and II, the two autistic children IV2 and IV4 shared a rare variant (ATP2B4 c.G1996A).

The 21 and 22 preliminary screening variants in autistic children IV2 and IV4, respectively, were verified by Sanger sequencing. The sequencing maps of three variants in IV2 and four variants in IV4 are provided in Fig. 2. As shown in Tables III and IV, these variants were ultimately identified as highly credible autism-related variants. The sequencing maps of the remaining variants will be provided upon reasonable request (see Availability of data and materials below). To determine the source of these preliminary screening variants, the sequencing data of the parents were analyzed. The results indicated that all the preliminary screening variants in the two autistic children were inherited from their parents, who were heterozygous for all of these variants, and there were no de novo variants (Tables III and IV). The sequencing maps of the parents will be provided upon reasonable request (see Availability of data and materials below).

A total of 11 candidate autism-related genes (indicated with ^a in Table III) in autistic child IV2 were screened according to the evidence score of >4 in AutismKB or ≥3 in AutDB. The candidate genes with an evidence score of >16 in AutismKB were regarded as highly credible autism-related genes, which included LAMC3, JMJD1C and CACNA1H (indicated in bold in Table III).

Based on the possible autism-related gene functions in Gene Ontology (GO; updated on January 1, 2018) database (https://www.uniprot.org/help/gene_ontology), LAMC3 was associated with astrocyte development, JMJD1C was related to transcriptional regulation and histone demethylase activity (that is, chromatin remodeling) and CACNA1H was involved in neuronal action potential and membrane potential regulation (that is, synaptic function).

The specific missense variant loci of these three genes were reported for the first time in the present study, to the best of our knowledge. None of the loci were found in AutismKB 2.0 autism database (version 2018), AutDB autism database (version 2020) or HGMD database (version 2019). The three gene variants were located outside the repeat region. The MAFs of LAMC3 c.G2687A were 0.001, 0.01 and 0.008 in East Asian populations (1000G, ExAC and gnomAD, respectively). However, CACNA1H c.A1046C and JMJD1C c.A4643G were not detected in these populations. According to the HGMD database, 40.91% (9/22), 92.00% (23/25) and 92.73% (51/55) of the missense variants in LAMC3, JMJD1C and CACNA1H genes, respectively, were pathogenic. The functional effects of the variants in these three genes, as predicted by 10 software packages are listed in Table V. Notably, eight of the 10 packages predicted LAMC3 c.G2687A as harmful (i.e. very high harmfulness), eight of the 10 packages are predicted CACNA1H c.A1046C as harmful (i.e. very high harmfulness) and two of 10 software predicted JMJD1C c.A4643G as harmful (i.e. high harmfulness).

In addition to LAMC3 and JMJD1C inherited from the father III3 (descended from I1 and I2), there was also the cumulative effect of CACNA1H inherited from the mother III4 (not descended from I1 and I2), i.e. there was a cumulative effect of CACNA1H with LAMC3 and JMJD1C. The results of the present study could provide guidance for the prenatal diagnosis of ASD. If the father III3 and mother III4 of autistic child IV2 give birth again, prenatal genetic diagnosis should be conducted on these three highly credible autism-related genes.

A total of 10 candidate autism-related genes (indicated with ^a in Table IV) in autistic child IV4 were screened according to the evidence score of >4 in the AutismKB or at ≥3 in the AutDB. The candidate genes with an evidence score of >16 in the AutismKB were considered as highly credible autism-related genes, which included SCN1A, SETD5, CHD7 and KCNMA1 (indicated in bold in Table IV).

Based on the possible autism-related gene functions in GO database, SCN1A was related to axon initiation segment and neuronal action potential transmission (synaptic function), SETD5 was involved in transcriptional regulation and histone acetylation regulation (i.e. chromatin remodeling), CHD7 was associated with central nervous system development (i.e. synaptic function) and transcriptional positive regulation, and KCNMA1 was responsible for the regulation of membrane potential (i.e. synaptic function).

SCN1A c.G1499A has been reported to be associated with Dravet syndrome in the HGMD database (updated in 2019), whereas the missense variant loci of the remaining three genes were first discovered in the present study, to the best of our knowledge, and were not found in the AutismKB 2.0 database (version 2018), the AutDB database (version 2020) or the HGMD database (version 2019).

All variants in these four genes were located outside the repeat region. The MAFs of SCN1A c.G1499A, SETD5 c.G3325A and KCNMA1 c.A34G were (0.003, 0 and 0.0016), (0.001, 0 and 0.0032) and (0.001, 0 and 0.0014) in East Asian populations (1000G, ExAC and gnomAD), respectively; whereas CHD7 c.G1892A was not found in these populations. According to the HGMD database, 51.62% (860/1666), 18.75% (6/32), 20.14% (171/849) and 66.67% (6/9) of the missense variants in SCN1A, SETD5, CHD7 and KCNMA1 genes, respectively, were pathogenic. The functional effects of the variants in these four genes predicted by 10 software are listed in Table V. Notably, seven of the 10 software packages predicted that SCN1A c.G1499A was harmful (i.e. very high harmfulness), nine of the 10 software packages predicted that SETD5 c.G3325A was harmful (i.e. very high harmfulness), seven of the 10 software packages predicted that CHD7 c.G1892A was harmful (i.e. very high harmfulness) and two of the 10 software packages predicted that KCNMA1 c.A34G was harmful (i.e. high harmfulness).

In addition to SCN1A, SETD5 and KCNMA1 inherited from the mother III6 (descended from I1 and I2), there was also the cumulative effect of CHD7 inherited from the father III5 (not descended from I1 and I2), i.e. there was a cumulative effect of CHD7 with SCN1A, SETD5 and KCNMA1. The results of the present study may provide guidance for the prenatal diagnosis of ASD. If the father III5 and the mother III6 of child IV4 give birth again, prenatal genetic diagnosis should be performed on these four highly credible autism-related genes.

The filtering procedure and number of the remaining variants in IV2 and IV4 autistic children are summarized in Table VI.