Hierarchical Genomic Analysis of Susceptibility to Schizophrenia in Sudanese multi-case families

Study design

This is a case pseudo-control, population-based study encompassing the genomic data of seven families diagnosed with schizophrenia by specialized psychiatrists based on DSM-IV and ICD-10 criteria, each with more than one affected member. Patients were selected for genetic analysis based on being siblings from the same family, with no age or gender specifications. Samples were collected from various psychiatric clinics and hospitals in Khartoum, Sudan, during the period from 2019 to 2021. The total number of participants is 29, including 18 cases and 11 controls.

Next generation sequencing

Whole Exome Sequencing (WES) was performed for one family, while Whole Genome Sequencing (WGS) was conducted for the other six families. The sequencing was done by BGI Tech Solutions (Hong Kong) Co., Limited. The sequencing service included generating at least 90 GB of data per sample, followed by data cleaning.

Basic Bioinformatics analysis was done by BGI analysis team encompassed several steps: data filtering, alignment of reads to the human reference genome (UCSC build HG19), then SNP and InDel calling and annotation of variants to corresponding gene functional units. Data received in Variant Call Format (VCF) files including information about Genes, Transcripts, Function, Variants type, Frequencies and others.

Genes filtering and variants selection

From the Variant Call Format (VCF) files of different variant types, including SNPs and INDELs, filtration of variants shared between all samples from all families, including patients and sibling controls, was performed using Linux. A family-based analysis was then conducted to compare variants shared between patients from the same family but not found in control siblings, using the R program. This comparison was based on variant transcripts or positions and was applied to all variants.

These shared variants, whether between all samples or only between patients, were further filtered according to their impact effect. The impact effects were classified as high, moderate, modifier, and low, based mainly on the variants’ gene functions, excluding intronic, intergenic, or pseudogene regions. Only regions and genes with high and moderate impact effects were included in the final analysis. SNPs and INDELs annotated and send from BGI Company using multiple annotation software for VCF files including clinvar, 1000G, cosmic and other databases.

In silico analysis

Five databases were used to analyze the function, clinical impact, and population allele frequencies of mutations: The National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/), Ensembl genome browser 108 (https://www.ensembl.org/index.html), the Genome Aggregation Database (gnomAD) (https://gnomad.broadinstitute.org/), Varsome (https://varsome.com/), and ACMG classification by The American College of Medical Genetics and Genomics (https://wintervar.wglab.org/).

Most of the genes filtered from the previous steps were checked for disease association. Genes already known to be associated with schizophrenia were identified using ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/), GeneCards (https://www.genecards.org/), the disease gene database DisGeNET (https://www.disgenet.org/), which contains data about human diseases from literature and GWAS catalogues, and Enrichr - Ma’ayan Laboratory, Computational Systems Biology (https://maayanlab.cloud/Enrichr/).

Each gene was then assessed for related biological pathways to identify any functions linking them to the pathogenicity of the disease using the Reactome database (https://reactome.org/).

The final set of genes was used to draw networks, including genes and protein interactions, using Cytoscape (https://cytoscape.org/), which relies on multiple network databases including GeneMANIA, Uniprot, Reactome, and BioGRID, as well as another network program, ConsensusPathDB (http://cpdb.molgen.mpg.de/).

Tissue expression for single or multiple genes was analyzed using GTExPortal (https://gtexportal.org/home/), providing details about tissue expression of these genes and detecting bulk tissue expression for specific genes accurately by sorting tissues based on median TPM (Median Transcripts Per Million) and displaying gene expression in a heat map view. The analysis is summarized in a flowchart (Figure 1).

Figure 1. Analysis flow chart.

Families

Form these families patients already diagnosed and on treatments, the main age group is (20-40) years old, there are no significant gender differences 10 females and 8 males.

Family 1 (F1): Both patients are males undergoing treatment. Their age at diagnosis was 17 and 18 years old. Their elder brother, who is healthy and free from mental disorders, served as the control. Additionally, two sisters from the same family have learning disabilities, with one showing typical hallucination and delusion symptoms but remained undiagnosed due to financial constraints. Both patients in F1 had poorly controlled symptoms despite medication.

Family 2 (F2): This family includes four affected members: two sisters, one brother, and their mother. Sequencing was successful for two affected sisters and two controls: their elder brother and younger sister. The two affected sisters’ age at diagnosis was 10 and 25 years old.

Family 3 (F3): Two affected members—a male and a female—and their healthy elder brother were sequenced.Their age at diagnosis was 35 and 25 years old, respectively. Both patients had poorly controlled symptoms despite antipsychotic medication.

Family 4 (F4): This family has two affected male patients. Their age at diagnosis was 20 and 24 years old, the elder patient also suffers from diabetes with associated complications, while the younger patient showed good compliance with medications. They have a control brother.

Family 5 (F5): Four affected siblings—three males and one female—and one unaffected sister as control. Their age at diagnosis was 17, 21, 23 and 27 years old. The father occasionally experiences delusions and paranoid ideas but has never been diagnosed with schizophrenia.

Family 6 (F6): Four affected female patients and their elder brother as control. Their age at diagnosis was 18, 33, 35 and 36 years old. The eldest sister showed severe symptoms with a long duration of illness and poor medication control, whereas the other patients had shorter duration and milder symptoms upon diagnosis.

Family 7 (F7): This family consists of two affected sisters as patients and their elder brother as control. They are the youngest siblings in the family, and both sisters have shown a good response to medications. Their age at diagnosis was 21, and 23 years old.

Classification of variants

Results were categorized and analyzed hierarchically based on their frequency within the families and their potential contribution to the phenotype (Figure 2). The flowchart summarizes data ranging from common variants shared between patients and family controls to rare variants.

Figure 2. Hierarchy representing the result ascending sequence trying to understand the phenotype.

The presentation of results depends on the degree of sharing between cases, with a focus on variants classified as having high impact effects. These variants are further categorized based on allele frequency, whether they are either common, or rare variants, with minor allele frequency (MAF) less than 0.01 or more than 0.05 (0.01< MAF <0.05) to assess their potential contribution to the phenotype.

Initially, all 29 samples, including both schizophrenic patients and sibling controls, were analyzed together. Assuming that the control sibling may carry the mutant allele without manifesting the condition. From the variants shared across all samples, 47 high impact SNPs and 42 high-moderate impact INDELs were selected. Most of these prioritized common variants are benign with some considered as likely pathogenic or of uncertain significance according to ACMG classification. Some already known to be associated with schizophrenia.

Analysis of the first level (common variants)

The analysis of shared SNPs among all samples, including both patients and sibling controls, identified 47 variants with high impact mutations. These mutations are associated with functions such as structural interaction variants, stop lost, splice donor, splice acceptor, and stop gained, suggesting their relevance to schizophrenia within these families.

Most of these variants exhibit a high frequency of the ancestral allele (1-0.8), with two exceptions noted in SND1 and NCAM1, where their frequencies are 0.6. NCAM1 is known to be associated with schizophrenia and is overexpressed in the frontal cortex of the brain.

The majority of these genes show expression in the brain, primarily in the cerebellum, frontal cortex, and spinal cord (Figure 3).

Figure 3. Major genes including (SNPs, INDELs) with different brain expression show cerebellum higher expression.

Symbols/Abbreviation: SNPs: Single nucleotide polymorphisms;

INDELs: Insertion/Deletion.

25 variants from the common variants group were identified despite being shared in all cases and the small sample size, their selection based on their high impact classification according to our sequencing annotation and their brain expression.

Most of these variants are found in genes show high brain expression, the majority expressed in the brain frontal cortex, while another eight genes are highly expressed in the cerebellum (Table 1).

Among these genes, DPYSL2, MPP2, ALDH3B1, NCAM1, CHRNA3 and GPRIN2 are known to be associated with schizophrenia, whereas NCOR2 is associated with bipolar disorder according to disease gene databases.

Comparison of the shared variant frequencies between the general population and Sudanese population showed almost no differences confirming that these variants are common in both populations. However, there are exceptions among four INDEL variants (rs11284059, rs11448549, rs34377180, and rs34373121), which show significant lower frequencies in the Sudanese genome (0.38-0.44) compared to the general population.

Family based analysis including rare variants: Variants collected from all families included 388 SNPs and 188 INDELs, and we selected 41 genes containing both SNP and INDEL variants. The selection was based on variants being shared between patients and absent in controls within each family, and then being shared between at least two families, regardless of having the same variant ID. Importantly, there were no variants shared between all seven families. There was a significant difference in the presence of high impact variants between cases and sibling controls in families, with a P value of 0.00862. From this group, we focused only on high-impact rare variants and further narrowed down to ultra-rare variants with an allele frequency ≤ 0.05%.

One INDEL variants shared between all samples including cases and controls presented in Table 1, irrespective of allele frequency, in NCOR2 (rs143952466) classified as rare variants with allele frequency assessed by gnomAD database (AF = 0.00002060) mainly expressed in Pituitary gland and Cerebellum. Associated with Bipolar Disorder according to diseases gene database. Classified by ACMG as uncertain significance.

Most of these rare variants are present in two or three families and are classified as likely benign/benign according to the ACMG classification. However, 12 variants were considered of uncertain significance and are located in genes associated with schizophrenia according to the OMIM and ClinVar 2019 databases. This group of genes is associated with Neurotransmitter Clearance according to the Reactome database, with a P value of 0.03475, suggesting a potential link to schizophrenia. However, the most highly involved biological pathway is the Biosynthesis of Maresins, which is more likely associated with type 2 diabetes.⁸

Other genes shared between families with more than two variants regardless of being rare variants, including both SNPs and INDELs, are CYP2C8, MUC6, PABPC3, COMT, and PADI4. Most of these variants are classified as benign by ACMG classification.

Genes known to be associated with schizophrenia include the COMT gene, which presents a shared mutation in four variants, one of which (rs199710929) is ultra-rare with an allele frequency of 0.0001402. The PADI4 gene has three variants (rs874881, rs2240335, rs2240335) with an allele frequency of ≤ 0.5%, and the CTNNA2 gene shows the same variants in two families (rs7609047, AF = 0.2435). The NAALAD2 SNP (rs139259253) has a rare allele frequency of 0.00008843 and is classified by ACMG as of uncertain significance. Additionally, the SNTG1 gene (rs770648099, AF = 0.2098) is highly expressed in multiple brain regions, with the highest expression in the brain’s frontal cortex.

The MUC19 gene presented with five different variants, including two SNPs, two INDELs, and one deletion (SVs). The most interesting variant is a rare SNP (rs73100424) with an allele frequency of 0.0415, which is functionally classified as of uncertain significance.

The total rare variants, some of which were mentioned previously, include ultra-rare variants with an allele frequency of ≤0.1%. None of these variants are pathogenic, but 12 are considered of uncertain significance according to ACMG classification. These variants are located in the following genes: ST6GALNAC1, COMT, CYP2D6, NAALAD2, CDNF, PRSS3, PRAMEF2, ZDHHC11, TAS2R46, CYP2C8, NBPF9 and OR5K4 (Supplementary Table).

Network analysis

The gene network includes shared genes between all samples with high impact common variants revealed 14 important central genes. Among these, SON and NCAM1 have been mentioned previously due to either their lower frequency or their likely pathogenic status. The main centrality genes include RIPK1, RBBP4, MAP 3K1, NCOR2, and ZNF2. Only SON, CHRNA3, and DPYSL2 are known to be associated with schizophrenia. Other genes in the network include TAAR9, CREB3L1, NET1, LTF, and VEGFC. Most of these variants are expressed in the brain, but unlike the known cerebral regions associated with schizophrenia, there is high expression in the cerebellum (Figure 4).

Figure 4. High impact variants associated network for SNPs and INDELs.

By using the STRING database and creating a more extensive network involving the previously mentioned genes, the expanded analysis revealed the involvement of additional genes. Together, these genes are related to three schizophrenia-associated biological pathways: catecholamine biosynthesis, metabolism of serotonin, and enzymatic degradation of dopamine by monoamine oxidase.

Other genes were studied because they are located at the center of each family network, and also appear in the network done by cystoscope joining all rare family variants. Genes were selected based on being shared between patients but not family controls. Three of these centrality genes are known to be associated with schizophrenia: PTPN6, which is central in two family networks with similar variants (rs2110072, AF = 0.9477) related to the immune system and cell communication; CTNNA2 (rs7609047, AF = 0.2534), associated with apoptotic cleavage of cell adhesion proteins, potentially contributing to neurodegeneration and neurodevelopmental issues; and GRIA2 (rs34460606, AF = 0.004842), which is central in the gene network and can contribute to the pathogenicity of schizophrenia due to its involvement in the activation of NMDA receptors, postsynaptic events, and neurotransmitter signaling pathways. Other genes from different family networks known to be associated with schizophrenia include HLA-B, HDAC9, PTPN1, TBP, and DGKH. Additionally, MAP 2K3 (rs55796947) and TMEM216 (rs10897158) both presented with SNPs that are not rare but are classified as likely pathogenic as germline mutation by Varsome data base.

In each family shared genes networks also in all families rare variants network, there is a cluster of ZNF genes interacting within a large group of the same gene family, but not interacting with other gene families. For the general rare variants network the ZNF717 was the centrality gene which is not known to be associated with schizophrenia.

Ethical approval and consent to participate

The project has been approved by the National Research Ethics Review Committee of the Health Research Council, Federal Ministry of Health, Republic of Sudan in 22/10/2019 under proposal no. (5-10-19) entitles (Whole genome sequencing in analysis of schizophrenia patients and their antipsychotic drug response variation in Sudan) submitted by Dr. Aza Elshafie Khidir Saeed, of the Institute of Endemic Diseases/University of Khartoum.

Patients were sampled and interviewed by the principal investigator following psychiatrist approval. Written informed consent was obtained from the patients themselves or their guardians after providing an informative description and discussion and approved by National Research Ethics review Committee of the Health Research Council, Federal Ministry of Health, Republic of Sudan.