Population-based study of genetic variation in individuals with autism spectrum disorders from Croatia

Autism Spectrum Disorders (ASD, MIM209850) are a severe neuropsychiatric disorders, primarily characterized by abnormalities in social behavior, communication and language, with patterns of restricted and repetitive interests [1]. Commonly associated symptoms include aggressive and self-injurious behaviors, anxiety, heightened sensitivity to stimuli, and seizures [2, 3]. Severity of symptoms in ASD can vary widely among cases, from individuals with mental retardation and no language production to relatively high-functioning individuals with normal to superior intelligence but significant difficulties with social interaction [4–6].

Multiple lines of evidence converge to suggest that ASDs represent the most heritable neurodevelopmental and psychiatric conditions. Evidence from twin and family studies suggests that the rate of autism in siblings of affected individuals is 2-6% [7], with 92% concordance of ASD in monozygotic twins and over 10% concordance in dizygotic twins [8]. Previous linkage and candidate gene studies have identified several chromosomal regions with autism susceptibility loci [3, 9–13]. A recent study identified a common genetic risk factor underlying ASD employing a genome-wide association strategy in 780 families (3,101 subjects) with affected children, and a second cohort of 1,204 affected subjects and 6,491 control subjects; these findings were replicated in an independent study cohort of 483 ASD families [14]. Furthermore, substantial progress towards the identification of genetic risk variants has come from recent characterization of structural variation (i.e., copy number variation or CNV). Studies of CNVs in a cohort of 859 ASD cases and 1,400 healthy control children revealed structural variants in previously reported ASD candidate genes (NRXN1 and CNTNAP4) and in multiple novel susceptibility genes encoding cell-adhesion molecules, including NLGN1 and ASTN2. Furthermore, CNVs observed in cases but not controls were found within and surroundings genes involved in the ubiquitin pathway (UBE3A, PARK2, RFWD2 and FBXO40) [15]. In a related study, we performed high-density genotyping and obtained data on 3,832 individuals from 912 multiplex ASD families in the Autism Genetic Research Exchange (AGRE) and 1,489 unrelated control subjects. Through prioritization of exonic deletions (eDels), we recovered genes in which structural variants were present in multiple unrelated probands but not unrelated controls. The 42 genes identified by this method include NRXN1, UBE3A, MADCAM1 and BZRAP1 [16]. Studies of large families with shared ancestry have reported several other autism loci including large, inherited, homozygous deletions in neuronal cell-adhesion genes, for example PCDH10 (protocadherin 10) [17]. These findings highlight the utility of "homozygosity mapping" and a need to search for additional large families across diverse populations.

Based on these breakthroughs and lessons learned from studies of other complex diseases (i.e. diabetes and breast cancer), it is becoming apparent that there is a need to assess large sample collections that number in the tens of thousands. Population structure poses a challenge for large-scale disease-association studies; mismatched ancestry between cases and controls in a genome-wide association studies is a potential source of spurious associations [17]. These needs and challenges bring together usually separated fields of genetics, epidemiology and studies of human diversity. Also, it is important to relate genetic findings found in patients from developed countries to the studies of diverse populations in less developed parts of the world. Currently DNA has been assembled from only a few thousand individuals with autism worldwide. Significant advances in genetic analyses of ASD in diverse populations will require the collection of samples across many populations using rapid and affordable screening and diagnostic tools.

Isolated populations have contributed to the discovery of genetic factors for many Mendelian disorders [18]. It has been argued that isolated populations may facilitate the identification of susceptibility loci for complex diseases because of a reduced genetic diversity, i.e. reduced genetic heterogeneity [19, 20]. Although several Croatian island populations have been described as genetic isolates, the Croatian population as a whole does not represent a population isolate as conventionally defined [21, 22]. However, Croatia is a small country with a low migration rate, in which large extended families often live in the same household or same neighbourhood and, therefore, is perfectly suited for family-based genetic studies. We expect that a subset of affected individuals may segregate a limited number of causative ASD susceptibility alleles (causal variants); in some cases these loci or disease-associated haplotypes identified in genetic isolates may differ from those in more genetically diverse populations, such as Europeans in North America. We also expect that a subset of rare variants associated with ASD susceptibility (SNPs and CNVs), which were identified in a large set of Europeans from North America, may be more prevalent in specific ethnic groups like the Croats. Therefore, to improve our understanding of global patterns of human genetic variation and set the stage for large-scale disease studies, efforts are underway to evaluate genetic variation using high-density genotype data for human populations worldwide [23–25].

With the intention to use the Croatian population in future genome-wide association studies, we performed a genome-wide high-density genotype pilot study of 103 ASD cases and 203 control samples to characterize the genetic structure of the population. Our analysis focuses on (1) the genetic diversity of the Croatian population, and its relationship to other populations, (2) the degree of relatedness of the population through the analysis of Runs of Homozygosity (ROHs), and (3) the distribution of large (>500 Kb) copy number variations, a major mechanism for structural variation in human genomes, and its possible association with ASD.

To understand the genetic homogeneity of the Croatian cohort and its position in the worldwide human population, we used phylogenetic tree analysis, principal component analysis (PCA) (not shown), and multi-dimensional scaling (MDS), three tools commonly used in the study of human population genetic diversity [23, 24]. We used data from the Human Genome Diversity Project (HGDP), which consists of genome-wide SNP data of 1007 individuals from 51 populations, divided into 7 geographic regions. Of these individuals, 153 were from 8 populations from the European region. We plotted the results of MDS and tree analyses using control subjects from Croatia and the HGDP samples (Figure 1). All three methods show the same overall pattern emerging from previous HGDP studies: the relative positions of the major regions largely correlate with their relative geographical locations, and reflect the human migration history well. Humans migrated from Africa to the Middle East (Figure 1A-B); the first branching took place, and one group migrated to Europe, while the other migrated to South Asia. A second major branching then occurred, with one migration to America (through Alaska), the other countinuing to East Asia and eventually Oceania. As expected, the Croatia control samples are entirely contained within the European samples. A focused analysis of Europeans in the HGDP panel (Figure 1C-D) with the Croatia control samples situated in the middle of phylogeny or MDS plots shows the relationship of this subpopulation to other European populations. The Sardinians and the Basques seem to be the farthest, with Italian samples placed in between them and the Croatian controls; the other subpopulation that is separate from the Croatia controls is the Adygei. The Croatians are most similar to the Russians and the Orcadians. Both the PCA (not shown) and MDS plots show the relative diversity of the Croatian controls (indicated by the area spanned by the samples on the MDS and PCA plots) to be similar to the other populations. The boxplot of pairwise Identity by State (IBS) distances for European populations (Figure 2A) showed the Croatian control cohort has divergence similar to other European subpopulations (except Basques and Sardinians who are genetically more homogeneous), and smaller than the European population as a whole.

We performed MDS on the Croatian controls and cases, as well as NINDS (National Institute of Neurological Disorders and Stroke) control cohorts to check the relationships between the three cohorts (Figure 2B), and found that the Croatian control and cases are well mixed but separate from the NINDS controls. Randomization tests by shuffling between Croatian control and NINDS controls showed the median pairwise IBS (identity-by-state) distance across the two subsets (0.28) is significantly higher than expected by chance (p = 0.00273, see Figure 3); this agrees with the MDS visualization that the two cohorts are genetically separate. Analysis of fixation index (F_ST) led to the same conclusion (see Table 2 for top regions with large F_ST). In summary, the Croatian and the NINDS cohorts have significant genetic difference, and care must be taken to correct for population stratification in combining the two datasets for association studies.

Motivated by the study demonstrating the positive correlation between runs of homozygosity and age in a North American cohort of European descent due to urbanization [35], we examined age correlation with ROH frequency and size in the 103 ASD subjects (Figure 4 and Table 4). We found only borderline significance between age and average ROH size, though the lack of significance may be due to the fact that the Croatia ASD cohort is 40 years younger (mean age is 21.5 y as opposed to 61.7 y in Nalls et al. [35]).

ROH regions were common or specific for a demographic group (Croatia vs. North America) or a disease status (ASD vs. controls). Although the longest ROH regions do not need to necessarily coincide with the disease susceptibility, we cannot exclude a possibility that recessive variants, such as single nucleotide or copy number variants within these regions contribute to disease susceptibility. To identify which regions contribute to the significant difference between Croatia case and control groups, we sorted the 7,740 endpoints of all autosomal ROHs across Croatia case and control groups, computed the significance of association between ROH status and case/control status by Fisher's test for each pair of adjacent ROH endpoints, and examined the significance of each such ROH interval where at least 5% (15) individuals have the ROH (Table 5). The most significant ROHs are in two chromosomal regions: chr2:82,040,863-82,677,336 (2p12) and chr6:29,481,192-29,719,410 (6p21.1) (Figure 5 and Figure 6). It is striking these two homozygosity regions are near several genes previously associated with neurological and psychiatric disorders. The chromosome 2 interval (636 Kb) is located downstream of CTNNA2 (alpha catenin) and LRRTM1 (leucine-rich repeat transmembrane neuronal protein 1), and contains a segment (between SNPs rs2685729 and rs2862972) significantly enriched in ASD subjects (top p = 2.5 × 10^-3). The LRRTM1 gene is associated with left-handedness (patients with autism and many other psychiatric disorders have a higher tendency to be left-handed [39]) and may increase risk of schizophrenia [40, 41]. The chromosome 6 interval (238 Kb) is a gene-rich region located distal to the major histocompatibility complex. This ROH is shared by 17 ASD subjects (16.5%) and 19 controls (9.6%) spans 8 genes, including a cluster of olfactory receptors (OR11A1, OR10C1, OR2H1, OR2H2), the MOG locus, and the GABBR1 locus. The MOG gene (myelin oligodendrocyte glycoprotein) is a transmembrane protein expressed on the surface of myelin sheath, and may be involved in multiple sclerosis (OMIM#126200) [42, 43]; recent studies also found that MOG is differentially expressed in prefrontal cortex and temporal lobes in schizophrenic patients [43, 44]. The most significant interval (between SNPs rs1233399 and rs29225; top p = 3.4 × 10^-3) spans the 3' end of the GABBR1 (Gamma-aminobutyric acid (GABA) B receptor, 1) and the overlapping UBD (Ubiquitin D) loci. GABBR1 is highly expressed in nervous systems and has been suggested as a candidate region for neuropsychiatric disorders [45, 46]. A recent study showed that the expression levels of GABBR1 are significantly lower in several brain regions in patients with autism compared with healthy controls [47]. Finally, we examined the median IBS distances between the subjects with regions of homozygosity on chromosomes 2 (13 Croatian cases and 10 Croatian controls) and chromosome 6 (17 Croatian cases and 19 Croatian controls) with permutation tests, but found that these subjects were not more closely related than the population average. Thus, the increased homozygosity in Croatia ASD cases cannot be explained by inbreeding among close relatives in a small number of selected families. However, the observed homozygosity may be the result of shared ancestral alleles from generations back in time.

CHR	START	END	SIZE (KB)	#SNPs	#ROH&Case	#Case	#ROH&Ctrl	#Ctrl	#ROH	#Samples	OR	P
2	82040863	82067472	26610	3	11	103	8	198	19	301	2.83	0.0424
2	82067472	82143271	75800	11	11	103	7	198	18	301	3.25	0.0195
2	82143271	82146419	3149	3	12	103	7	198	19	301	3.58	0.0105
2	82146419	82235100	88682	18	12	103	8	198	20	301	3.12	0.0153
2	82235100	82276290	41191	6	12	103	9	198	21	301	2.76	0.0305
2	82276290	82330384	54095	4	12	103	8	198	20	301	3.12	0.0153
2	82330384	82342863	12480	2	13	103	8	198	21	301	3.42	0.0081
2	82342863	82367354	24492	7	12	103	8	198	20	301	3.12	0.0153
2	82367354	82677336	309983	20	12	103	7	198	19	301	3.58	0.0105
2	82677336	82706135	28800	3	12	103	5	198	17	301	5.06	0.0025
2	135475565	135478814	3250	2	4	103	21	198	25	301	0.34	0.0488
4	34049422	34100634	51213	8	30	103	37	198	67	301	1.78	0.0421
4	34100634	34156199	55566	7	30	103	36	198	66	301	1.85	0.0393
4	34156199	34179397	23199	5	29	103	35	198	64	301	1.82	0.0386
4	34378375	34397251	18877	3	26	103	30	198	56	301	1.89	0.0420
4	34531827	34538866	7040	2	17	103	16	198	33	301	2.24	0.0326
5	42508871	42588188	79318	9	11	103	8	198	19	301	2.83	0.0424
6	29481192	29574935	93744	17	16	103	15	198	31	301	2.24	0.0442
6	29574935	29591890	16956	7	16	103	14	198	30	301	2.41	0.0255
6	29591890	29602876	10987	3	15	103	13	198	28	301	2.42	0.0348
6	29602876	29607484	4609	2	15	103	11	198	26	301	2.89	0.0156
6	29607484	29623781	16298	3	16	103	11	198	27	301	3.11	0.0056
6	29623781	29638829	15049	5	16	103	10	198	26	301	3.44	0.0040
6	29638829	29639232	404	3	15	103	9	198	24	301	3.56	0.0034
6	29639232	29641274	2043	2	14	103	9	198	23	301	3.29	0.0101
6	29641274	29645954	4681	4	14	103	10	198	24	301	2.95	0.0131
6	29645954	29647461	1508	2	14	103	9	198	23	301	3.29	0.0101
6	29647461	29666169	18709	6	14	103	8	198	22	301	3.72	0.0042
6	29666169	29689020	22852	11	13	103	8	198	21	301	3.42	0.0081
6	29689020	29716673	27654	2	12	103	7	198	19	301	3.58	0.0105
6	29716673	29719410	2738	2	11	103	7	198	18	301	3.25	0.0195
6	29719410	29723801	4392	3	10	103	7	198	17	301	2.92	0.0356
6	61949597	62025463	75867	2	11	103	7	198	18	301	3.25	0.0195
6	62030461	62041958	11498	4	12	103	9	198	21	301	2.76	0.0305
7	69474084	69475255	1172	2	2	103	16	198	18	301	0.23	0.0392
8	51323785	51339189	15405	2	7	103	30	198	37	301	0.41	0.0416
12	33162179	33170263	8085	2	10	103	6	198	16	301	3.43	0.0266
12	33170263	33329348	159086	27	10	103	7	198	17	301	2.92	0.0356
12	33329348	33331881	2534	2	11	103	7	198	18	301	3.25	0.0195
12	33331881	33335472	3592	2	12	103	8	198	20	301	3.12	0.0153
12	33335472	33349548	14077	2	13	103	11	198	24	301	2.45	0.0426
12	49797511	49861882	64372	6	10	103	7	198	17	301	2.92	0.0356

Three of the ASD cases have multiple large rearrangements; one adult male harbors two chromosomal rearrangements (separated by 117 Mb) on chromosome 1, a 335 Kb deletion on 1p35.3 and a 882 Kb duplication on 1q21.1 (Figure 7a). In a second individual (adult female) we identified two distinct large CNVs: a 838 Kb deletion in an intergenic region between GRM8 and POT1 on 7q31.33 (chr7:124,752,465-125,591,197) and a 691 Kb deletion at 22q13.33 (chr22:48,833,840-49,524,956). This deletion encompasses the SHANK3 gene and cytogenetic rearrangements affecting this gene have been previously associated with ASD [60, 61]. It has been previously reported that the 22q13 deletion syndrome, also known as Phelan-McDermid syndrome (MIM606232), is characterized by severe neonatal hypotonia and global developmental delay, normal to accelerated growth, absent to severely delayed speech, and minor dysmorphic features [62, 63]. The third individual (male child) with multiple large rearrangements harbors a 534 Kb intergenic deletion at 2q21.1 (chr2: 129,716,352-130,250,751) and a 400 Kb duplication at 3q29 (chr3:195,987,870-196,387,903). These large structural variants, especially those found on multiple chromosomes illustrate the need to combine high-density genotype analysis with karyotype analysis to test for the presence of balanced translocations and inversions.

In contrast to large rearrangements encompassing >10 Mbs that were observed in ASD cases, the largest rearrangement observed in control subjects is a 2.5 Mb duplication at 15q24.1-q24.2 (chr15:70,783,089-73,316,235), a known microdeletion/duplication region [64]. Also, we found one control subject to have two duplications within the DiGeorge (MIM188400) and Velocardiofacial (MIM192430) syndrome region on chromosome 22 (22q11.21), a known pathogenic rearrangement hotspot where large (>1.5 Mb) deletions are pathogenic [65]. These duplications (922 Kb in chr22:17,690,812-18,613,429 and 563 Kb in chr22:19,102,598-19,665,559) are separated by 489 Kb. The observation that we see more large CNVs associated with ASD subjects (>2 Mb: 5 cases, 1 control; Fisher's test; p = 0.019) suggests these large rearrangements may be rare ASD susceptible variants; follow-up investigations on larger cohorts will be necessary.

In this paper, we report analysis of a genome-wide high-density genotype study for ASD in Croatia. The dataset, although small, provides a distinctive opportunity to evaluate structural variation associated with a complex disease such as ASD in a genetically distinct population. In addition to the analysis of ASD, high-density genotype data from a well-defined geographic area of Croatia permitted also the study of the subjects' genetic substructure and diversity. We asked whether the subjects included in this study which were collected throughout Croatia, show a lower level of genomic diversity in comparison to subjects from other defined demographic areas. Addressing this issue may suggest that a set of 100 Croatian ASD subjects harbor a smaller set of susceptibility loci in comparison a randomly selected set of ASD subjects from a larger area. The Croatian samples represent a genetic continuum of relatedness with a subset of individuals closer to those from several regions in Europe, probably due to numerous historic migrations. The impact of fine-scale population history on the genetic substructure of individuals from Croatia has recently been evaluated using a small set of selected loci [21]. To our knowledge, this study is the first to provide a genome-wide view of the genetic composition of the Croatia population.

Although our study represents the first high-resolution genomic analysis of subjects with ASD in Croatia, several caveats need to be addressed. Apart from an obviously limited sample size, the classification of ASD subjects in this pilot project has been established using the DSM-IV criteria. Due to phenotypic heterogeneity of ASD and our finding of a wide-range of novel, and in some cases, exceptionally large CNVs, there is a need for accurate clinical diagnostic criteria for the diagnosis of ASD and related neurodevelopmental syndromes. Currently no other diagnostic tools were available in Croatian; in order to compare our findings to those in other large-scale genetic studies it will be necessary to translate gold standard diagnostic tools such as ADOS [66] and ADI-R [67] to Croatian and many other languages.

We used PennCNV to call copy number variations, though other CNV calling programs such as QuantiSNP [68] and cnvPartition (a built-in software tool in the Illumina BeadStudio software)are available An advantage of using multiple CNV calling programs is to reduce false positives, namely one can focus on CNVs returned by multiple programs and expect them to be "more reliable". Since our study focused on large CNVs (>500 Kb) we were able to validate these variants by visual inspection of the raw intensities(data not shown).

We found that runs of homozygosity are more frequent in Croatia ASD subjects than either Croatia or NINDS control cohorts. Given the low frequency of subjects with ROH, larger follow-up studies are necessary to confirm the findings. Our findings suggest there is a high chance that homozygosity mapping based on family studies can discover rare recessive variants. Our finding that several ASD subjects harbor exceptionally large chromosomal rearrangements previously associated with distinct clinical syndromes such as Velocardiofacial, Jacobsen, Phelan-McDermid, or 3q29 microdeletion/duplication syndromes (MIM609425; MIM611936) lends support for cytogenetic or high-density genotype testing for ASD subjects with complex clinical manifestation [69]. Furthermore, recent findings that seemingly diverse neurodevelopmental and neurological disorders are often associated with the same chromosomal rearrangement(s) provide support for in-depth phenotype analysis and longitudinal clinical evaluation of entire families. These family-based studies may also reveal whether observed structural variants are de novo events or rather inherited and present in additional family members. We suggest that deeper population screens combined with family-based genetic analyses will lead to improved understanding of joint actions of associated loci and molecular basis of neurodevelopmental disorders in populations worldwide.