Identifying Genomic Variations in Monozygotic Twins Discordant for Autism Spectrum Disorder Using Whole-Genome Sequencing.

Autism spectrum disorder (ASD) presents a set of childhood neurodevelopmental disorders with impairments in social communication and restricted, repetitive, and stereotyped patterns of behavior. Here, based on the whole-genome sequencing (WGS) data of three monozygotic twins discordant for ASD, we explored multiple patient-specific genetic variations and prioritized a list of ASD risk genes. Our results identified DVMT (discordant variation in monozygotic twin) observed in at least two twin pairs, including 14,310 SNPs, 2,425 indels, and 16,735 CNVs, referring to a total of 2,174 genes, and 37 of these were covered by all three types of variations. Gene ontology (GO) enrichment analysis of biological processes for 2,174 genes showed that the majority of these genes were related to neurodevelopmental processes. In addition, functional network analysis showed that there was a strong functional relevance between 37 genes covered by all three types of variations. In conclusion, for the first time, we conducted a comprehensive scan of genomic differences between monozygotic twins discordant for ASD, providing researchers with in-depth directions. It may also provide effective strategies for clinical treatment of individuals affected by ASD.

Introduction

Autism spectrum disorder (ASD) represents a set of childhood neurodevelopmental disorders with impairments in social communication and restricted, repetitive, and stereotyped patterns of behavior. It is a lifelong developmental condition with a great degree of heritability and heterogeneity. A recent nationwide population-based study estimated that ASD prevalence was at 2.47% among U.S. children and adolescents in 2014–2016. Particularly, the prevalence rate for boys is nearly five times higher than that for girls. Although ASD typically occurs in childhood, patients may even need nursing during their whole life, bringing a heavy burden to the whole society. ASD is now considered as one of the world’s most urgent public health issues. Since ASD is highly heritable, family studies reconfirm the strong genetic contributions to this disease, declaring it as one of the most heritable neurodevelopmental and neuropsychiatric disorders.2, 5 Growing evidence shows the importance of complex genetic factors in ASD (single-nucleotide variants [SNVs], indels [insertions or deletions], and CNVs),6, 7 and suggests the possible existence of hundreds of autism risk genes.8, 9, 10, 11 Whole-exome sequencing (WES)9, 12, 13, 14 and whole-genome sequencing (WGS)15, 16, 17 have been used to detect de novo mutations as long as inherited mutations were involved in ASD and other neuropsychiatric disorders, facilitating the discovery of rare pathogenic variations and risk genes. The curation of large genomic datasets has produced a collection of many genes implicated in ASD risk and ASD-linked genes.8, 10, 18 However, even considering all of the data from these complementary approaches, the effect remains incomplete, and up to now, only about 20% of the ASD cases have been explained.6, 19, 20

Monozygotic or “identical” twins have been widely studied to dissect the relative contributions of genetics and environment to human diseases, and they are widely used in traditional experimental research. While recent high-throughput research only concerned fraternal twins, few studies have empirically tested for genetic differences between monozygotic (MZ) twins discordant for ASDs. MZ twins were thought to originate from a single fertilized zygote and thought to be genetically identical. However, some MZ twins have the same genetic foundation but presented different phenotypes; in this case, MZ twins appear to have discordant phenotypes for ASD. To unearth certain causes of this difference, we conducted a comprehensive scan of genomic differences between MZ twins discordant for ASD based on WGS (Figure 1). In this study, we first collected three pairs of MZ twins discordant for autism and captured all types of genetic variation, including SNVs, small insertions or deletions (indels), and copy number variations (CNVs). Further, for each twin pair, variations that only presented in individuals with ASD were filtered as discordant variations in MZ twins (DVMTs), and DVMTs in at least two twin pairs were identified for further analysis. At the same time we got a gene set of 2,174 genes involved in all DVMTs, gene ontology (GO) enrichment analysis found many genes were related to the development of autism disease. In addition, functional network analysis showed that there was a strong functional relevance between 37 genes covered by all three types of variations, showing a preference for a diagnosis of autism.

Figure 1

The Bioinformatics Analysis Pipeline

The process includes three steps: the preprocessing of whole-genome sequencing data, the discovery of variants, and the integrative analysis of variant-related genes.

Results

Characteristics of WGS Data

After the filtering of raw data, the basic statistics collection of WGS data for three pairs of MZ discordant twins was concluded (Table S1). In summary, the average number of reads for raw data is 1,151,097,111. After stringent quality control at the sample level, the mean rate of clean data to raw data is 95.11%. Then we computed the genomes from 6 subjects and got an average of 99.53% coverage and 30.85× sequencing depth when mapping the clean data to hg38 with Burrows-Wheeler Aligner (BWA) (Table 1). Likewise, 99.25% of the genome was covered with at least 4× sequence depth. On average, the mapped reads for each sample were 1,040,389,962, and the mapping rate is 98.30%. The high genome coverage of WGS allows optimal variations calling and potentially captures many clinically relevant mutations.

Table 1

The Basic Statistics of Mapping and Quality Control

Twin Pairs	Total Reads	Mapped Reads	Mapping Rate (%)	Average depth (×)	Coverage ≥ 1× (%)	Coverage ≥ 4× (%)
Twin Pair 1
Twin 1a	1,051,026,934	1,034,220,852	98.40	30.63	99.42	99.15
Twin 2	1,079,153,824	1,062,219,892	98.43	31.34	99.40	99.11
Twin Pair 2
Twin 1a	1,042,955,285	1,025,555,250	98.33	30.35	99.41	99.12
Twin 2	1,094,044,728	1,070,511,132	97.85	31.87	99.43	99.17
Twin Pair 3
Twin 1a	1,074,131,972	1,057,584,856	98.46	31.38	99.75	99.48
Twin 2	1,009,031,528	992,247,792	98.34	29.57	99.76	99.49
Mean	1,058,390,712	1,040,389,962	98.30	30.85	99.53	99.25

The twin with autism in the twin pair.

Identifying DVMTs in Three MZ Discordant Twins

For SNVs and short indels (Figure 2), all SNVs and indels of each individual were first captured by the Genome Analysis Toolkit (GATK). At the same time, those variations in Single Nucleotide Polymorphism Database (dbSNP) were filtered out. Second, for each twin pair, variations that only presented in the individuals with ASD were identified, including 50,644, 42,387, and 49,869 SNVs and 11,118, 9,373, and 10,141 indels. For CNVs (Figure 2), DELLY was performed to identify duplication and deletion, respectively. CNVs with lengths in the range of 100 bp to 1 million bp were retained, including 50,644, 42,387, and 49,869 for duplication and 11,118, 9,373, and 10,141 for deletion for each twin pair. All these variations were observed only in individuals with ASD and were defined as DVMTs.

Figure 2

The Filtration Process of SNVs, Indels, and CNVs

Further, DVMTs that appeared in at least two twin pairs were identified (Tables S2 and S4). In particular, two CNVs with at least a 50% reciprocal size overlap with each other were the same, and this overlapping region was extracted as a CNVR (copy number variation region). Considering that, 14,310 CNVRs for duplication and 2,425 CNVRs for deletion were extracted (Table S6). Finally, those DVMTs were annotated using the ANNOVAR tool (Figures 3A–3D and 4A–4C; Tables S3, S5, and S7). We can see that variations located in the intergenic region account for the biggest proportion, the result of both SNVs and indels, in which non-coding RNA (ncRNA)_intronic has the biggest part of ncRNA. On the other side, no indels located in the exon and 5′ UTR compared to SNVs. For the final set of SNVs, we computed the mutation type distribution and found that C > T (G > A) accounted for the largest proportion (Figure 3E). Among those, transition has a higher rate than transversion (Figure 3F). For CNVs, we showed the distribution of duplication and deletion in three twin pairs, respectively (Figure 4B). Finally, to illustrate all kinds of variations, we displayed them in a Circos plot (Figure 5).

Figure 3

The Detection Results of SNVs and Indels

(A) Annotation result of SNVs. (B) Detailed distribution of ncRNA in SNVs. (C) Annotation result of indels. (D) Detailed distribution of ncRNA in indels. (E) The mutation signature of SNVs. (F) Distribution diagram of transversion and transition.

Figure 4

The Detection Results of CNVs

(A) Annotation result of DEL in CNVs. (B) The distribution of DUP and DEL in CNVs for each twin pair. (C) Annotation result of DUP in CNVs. DEL, deletion; DUP, duplication.

Figure 5

The Display of Circos for All Variations

From the innermost band to the outermost band, CNV (Deletion), CNV (Duplication), the density distribution of indels, and SNVs are indicated.

Prioritization of ASD Risk Genes Based on DVMT

While SNVs and indels usually cover one single gene, CNVs are typically large and may contain multiple genes, preventing the identification of the rightful pathogenic gene. To compile all possible candidate ASD risk genes and evaluate the influence of all variations on autism, we focused on a total of overlap 2,174 DVMT-related genes (Table S8), including 823 genes for SNVs, 557 genes for indels, and 1,174 genes for CNVRs (Figure 2). Considering the distribution of genes by ANNOVAR, 239 of these genes located on exons. Functional enrichment analysis was performed for these 2,174 candidate genes, showing that the majority of these genes were involved in the regulation of cell morphogenesis, the development of dendrite and synapse, the regulation of synapse, and modulation of chemical synaptic transmission (Figure 6A).

Figure 6

The Functional Analysis of Genes

(A) Enrichment of the biological process for 2,174 genes. (B) Network analysis of 37 genes by GeneMANIA.

To further explore candidate genes and prioritize important genes, we focused on 37 genes that appeared in all three types of variations (Table S9). We analyzed the functional network of these genes by GeneMANIA. Results indicated that 26 of 37 genes were contained in the network (Figure 6B), covering networks of co-expression, co-localization, and shared protein domain and genetic interactions, referring to 97.41%, 1.93%, 0.45%, and 0.21% of all genes, respectively, showing that there is a functional association among them. Referring to the network of 37 genes, several genes are related to important functions, which are also core connectors in the network. Among those genes, some are directly associated with ASD, such as Duchenne muscular dystrophy (DMD), a widely studied molecular factor that can be applied in the early diagnosis of autism; it contributes to the organization of autism-associated trans-synaptic neurexin-neuroligin complexes and to the clustering of synaptic GABAA receptors. Likewise, NRXN3, a type of neurexin gene that is a presynaptic cell-adhesion molecule affecting the function of synapses and mediating the conduction of nerve signals, plays an important role in normal brain development. ASMT might be a susceptibility gene for autism as indicated by a large-scale study of Chinese Han individuals, using target sequencing. ZNF713 and PTPRN2 are identified as new candidate susceptibility genes.30, 31 GPC5 suggests a possible role in growth regulation and in cognitive development.32, 33 Besides, those genes were related to a series of biological processes, such as growth regulation, cognitive development, normal brain development, and conduction of nerve signals.28, 29, 30, 31, 32, 33, 34 These functions demonstrate potential links in the progression of autism disorder disease.

Discussion

Growing evidence shows the importance of complex genetic factors in ASD (SNVs, indels, and CNVs),6, 7 suggesting the possible existence of hundreds of autism risk genes.8, 9, 10, 11 There exists several databases curating general autism risk genes, such as AutDB, autismKB, and NPdenovo. However, more studies are still needed to further refine and contribute to the existing genetic architecture of ASD,38, 39, 40 for the genetic contribution can only explain less than 20% of the patients. In this study, we reported the full detection of genomic differences between MZ twins discordant for autism and aimed to test whether somatic mutations play a role in autism etiology. Three genomic changes dominant in affected individuals were identified by deep sequencing in MZ discordant twin pairs. Our goal was the identification of mutations that were present only in the affected twin, which pointed to novel candidates for autism susceptibility loci.

Careful clinical characterization has revealed that ASD is extremely heterogeneous; variants may be linked to ASDs and other related disabilities, such as intellectual disability, or medical co-morbidities, such as seizures. Further investigation of clinical relevance to these variants is ongoing. To discuss the potential implications of our study, as a more general issue to the use of next-generation sequencing (NGS) for the study of autism and other neuropsychiatric disorders, we integrated evidence from protein-protein interaction (PPI) data curated from the String Database and collected other known neurodevelopmental disease risk genes from public databases (NPdenovo, AutDB, and autismKB) to find highly interactive genes, including a total of 4,434 genes, with an overlapping of 574 genes with our union gene set list (2,174 genes), and they are significantly enriched (p = 0.014, hypergeometric test). Based on PPI data, 20 of 37 genes are found to have a PPI with public risk genes. Evidence has shown that genes implicated in ASD also affect several neurodevelopmental disorders by disrupting a wide spectrum of pathways. For example, RBFOX1 plays a crucial role in the regulation of splicing and transcriptional networks in human neuronal development.

Genes displayed as core linkers in the network are also closely related to neurodevelopmental diseases (Figure 6B). After a systemic search of papers in PubMed, we found that several genes have been reported to be related to autism. Among those genes, DMD is widely studied in autism disease and shows a preference for a diagnosis of autism; it contributes to the organization of autism-associated trans-synaptic neurexin-neuroligin complexes and to the clustering of synaptic GABAA receptors. Particularly, males with DMD are at a very high risk of neuropsychiatric disturbance. Likewise, NRXN3 is a type of neurexin gene, and neurexins are presynaptic cell-adhesion molecules that affect the function of synapses and mediate the conduction of nerve signals and that play an important role in normal brain development and become candidate genes for autism. Wang et al. found that, among those, ASMT might be a susceptibility gene for autism, in a large-scale study of Chinese Han individuals, using target sequencing. ZNF713 and PTPRN2 have been identified as new candidate susceptibility genes.30, 31 GPC5 and DPP6 suggest a possible role of GPC5 in growth regulation and in cognitive development.32, 33

However, we have not yet proven whether these mutations, singly or in combination, contribute to the development of ASD. Given that this is only a pilot study for WGS in ASD and is based on a relatively small number of twins, the replication of our findings is needed to better describe the genetic architecture of ASD in the future. For other mutations, their role in the condition needs to be closely followed in the medical literature. The reported genetic variants also highlighted potential molecular targets for pharmacologic manipulation and open the way for personalized therapeutics in ASD. Promising results, thus far, indicate that WGS could eventually help direct personalized approaches to clinical management in individuals and families affected by ASD. In the future, results of WGS might allow earlier diagnosis of ASD, especially in siblings, for whom recurrence rates are approximately 18%.

Materials and Methods

Sample Selection

Three pairs of ASD-discordant MZ twins were recruited from the Child Development and Behavior Research Center (CDBRC), Harbin Medical University, Heilongjiang Province, China. More than two experienced psychiatrists independently issued ASD diagnoses according to the criteria of the Diagnostic and Statistical Manual of Mental Disorders, 5th edition (DSM-5). All probands were also administered by the Autism Diagnostic Observation Schedule—Generic (ADOS) and were all found to meet the criteria for autism. For each patient, we applied the ADOS module according to the level of expressive language. Cases with Rett syndrome, tuberous sclerosis, fragile-X syndrome, and any other neurological conditions suspected to be associated with autism were excluded by clinical examination and molecular genetic tests of the FMR1 gene. Approval from the Ethics Committee of Harbin Medical University (HMUIRB2012008) was obtained before initiating the study.

WGS

WGS was performed with HiSeq 2000 technology (Illumina) at the Beijing Genomics Institute (BGI) in Shenzhen.

Genomic DNA from 3 MZ twin pairs discordant for ASD was isolated from the peripheral blood and used to construct a genomic DNA library for WGS. DNA was extracted and sheared into fragments, which were then purified by gel electrophoresis. DNA fragments were ligated with adaptor oligonucleotides to form paired-end DNA libraries with an insert size of 500 bp. To enrich libraries for sequences with 5′ and 3′ adaptors, we used ligation-mediated PCR amplification with primers complimentary to the adaptor sequences. The DNA libraries were sequenced to generate 90-bp paired-end reads with at least 30× average genome coverage per sample.

Discovery and Characterization of Various Variants

Raw data were preprocessed by SOAPnuke to get clean data, which removed sequences of adaptor, low-quality reads, and unknown reads. Both raw data and clean data were quality controlled by FastQC. Filtered reads were then mapped to the human reference genome (hg38) with the BWA, v0.7.17. The BAM file was then sorted by genomic position and indexed using Sequence Alignment/Map tools (SAMtools). BAM files from the same sample were then merged into a single BAM file using MergeSamFiles in Picard tools. To prevent PCR artifacts from influencing the downstream analysis, we used MarkDuplicates in Picard tools to mark the duplicate reads, which were then ignored. Qualimap was used to calculate the basic statistics of BAM files.

Local realignment and quality recalibration were performed with the Genome Analysis Toolkit (GATK),v3.8. The HaplotypeCaller (HC) was used to detect SNVs and small indels (< 65 nt). We then used the GATK Variant Quality Score Recalibration (VQSR) to filter spurious SNVs and indels due to sequencing errors and mapping artifacts. We performed CNV analysis with DELLY. CNV calls with a p value < 0.05 and calls with >50% of q0 (zero mapping quality) reads within the CNV regions were removed (q0 filter). All variations are displayed in a Circos plot.

Variant Filtration and Risk Gene Prioritization

When assessing variants of SNVs and indels, we removed those variants common in the population (in dbSNP142) that were likely to be false positives. For each twin, we also removed variants of the normal sibling from the proband, and then we retained variants found in at least two pairs of twins. GATK SelectVariants was used to get intersect or union set from Variant Call Format (VCF) files. As for CNVs, sizes between 100 bp and 1,000,000 bp were preserved, since it is an interval region of genome; we considered CNVs with at least 50% reciprocal size overlap with each other to be the same and only retained the overlapped region for the following analysis, considered as the CNVR. We used BEDtools to achieve this goal.

The effects (e.g., missense, nonsense, or frameshift mutations) and classifications (e.g., in exonic, intronic, or intergenic regions) of variants across the genome were annotated by ANNOVAR. Finally, we got a gene set annotated from each type of variation and computed the union or intersect set of those genes.

Functional Analysis

For gene union set from three types of variation, we conducted a functional enrichment analysis of biological process by clusterProfiler. For gene intersect set from three types of variation, we explored the functional network of these genes by GeneMANIA.

Data and Software Availability

The sequencing data reported in this paper have been deposited at the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/; SRA: PRJNA471799).

Author Contributions

Y.H., Y.Z., Y.R., Y.Y., and J.Y. performed the analysis of sequencing data. Y.H., X.L., Z.G., and X.Z. collected the samples and performed the clinical examination and molecular genetic tests. Y.H., Y.Z., and Y.Y. wrote the initial manuscript. S.L. and D.W. improved the thesis of the section of the paper. L.W. designed the experiment. All authors reviewed the manuscript and provided comments that were incorporated in the final version of the manuscript.

Conflicts of Interest

The authors declare no competing interests.