A comparative bioinformatic analysis of C9orf72

Data collection

Genomic DNA, cDNA and protein sequences from a selection of species representing various taxonomic orders that have C9orf72 were selected for comparison from the Ensembl database (http://www.ensembl.org). Table 1 lists the species compared along with their gene and transcript IDs. Data describing the various features of the genes are collated in Table 2 to provide easy assessment of similarities and differences between the species.

Phylogenetic analyses

Ensembl Compara provides cross-species resources and analyses, at both the gene and amino acid sequence level. Ensembl Compara was used to create the phylogenetic tree from a multiple sequence alignment of protein sequences that are representative of the longest protein-coding translation of the gene from all species compared. The program TreeDyn (Chevenet et al., 2006); (http://www.phylogeny.fr) was used to visualize and annotate the phylogenetic tree.

Sequence alignments

The numbers and structure of alternative transcripts were acquired by searching the Ensembl database for each species of interest for C9orf72. Amino acid and cDNA sequences for C9orf72 from all the species studied were aligned using ClustalW2 (Larkin et al., 2007), a multiple sequence alignment software using default alignment parameters (http://www.ebi.ac.uk/Tools/msa/clustalw2/). BOXSHADE program was used to shade regions of similarity in the aligned sequences (http://www.ch.embnet.org/software/BOX_form.html). Matcher (one of the EMBOSS programs) was used to locally align shorter sequences (Rice, Longden & Bleasby, 2000; McWilliam et al., 2013).

Transcription factor binding site analysis

DiAlign TF and MatInspector (Cartharius et al., 2005) (http://www.genomatix.de/) were used to predict transcription factor binding sites within the nucleotide region encompassing the 5′-flanking region of the C9orf72 gene (2,000 bp upstream) up to the end of intron 1 of human, mouse and Fugu. DiAlign TF builds alignments from gap-free pairs of similar segments of the sequences, looks for local similarities and then searches for common regulatory sequences in those aligned regions. TF binding site matches were identified (in both aligned and non-aligned upstream region of the three species) by MatInspector using Matrix Family Library Version 11.0.

MEME (Bailey & Elkan, 1994), Tomtom (Gupta et al., 2007) and GOMo (Buske et al., 2010), three motif-based sequence analysis tools from the MEME suite of programs, were used to analyze the upstream regulatory region in human, mouse and Fugu. MEME (Bailey & Elkan, 1994) was used to identify new ungapped motifs of recurring fixed-length patterns in the chosen nucleotide region. The motifs identified were searched against different databases of known motif (transcription factor families) using Tomtom (Gupta et al., 2007). GOMo was used to identify GO (gene ontology) terms associated with the DNA regulatory motifs identified by MEME (Bailey & Elkan, 1994).

Comparative assessment of genomic synteny of C9orf72

To identify conserved syntenic regions, genes flanking C9orf72 in human, mouse and Fugu were annotated. Initial predictions were derived from the online server Cinteny (Sinha & Meller, 2007) (http://cinteny.cchmc.org/) and The Synteny Database (Catchen, Conery & Postlethwait, 2009). The programs identify syntenic regions across the genomes selected for comparison and analyse the extent of genome rearrangement using reversal distance as a measure. The predictions from these two programs were cross-referenced with those identified in the Ensembl database.

Gene location

The human C9orf72 gene is located on the short (p) arm of chromosome 9 (cytogenic position 21.2) open reading frame 72, from base pairs 27,546,542 to 27,573,863 (DeJesus-Hernandez et al., 2011; Renton et al., 2011). The gene, encoded on the minus strand of chromosome 9, consists of 11 exons that through alternative splicing are transcribed into 5 mRNA splice variants (DeJesus-Hernandez et al., 2011; Renton et al., 2011) (Fig. 1). The mouse gene, located on chromosome 4 on the reverse strand, like its human counterpart, also consists of 11 exons (Fig. 1). The Fugu C9orf72 gene differs from the human and murine gene in that it is encoded on the forward strand of the chromosome (scaffold_396) and comprises of 10 exons (Fig. 1). Although, Tetraodon has an even smaller genome than the Fugu, it was not considered here because the annotation of the Fugu gene entry in Ensembl is believed to be more complete.

Conservation of C9orf72 in vertebrates

Zhang et al. (2012) reported that C9orf72, a DENN-domain containing protein, could be traced all the way back to the last eukaryotic common ancestor. They also identified C9orf72 homologues in protozoans and other species outside of the metazoa. Since the focus of this report is on understanding the degree of conservation of C9orf72 amongst eukaryotes with particular emphasis on vertebrates, we looked at the phylogeny of C9orf72 from this perspective. The phylogenetic tree we constructed for the analysis of the evolution of C9orf72 gene amongst eukaryotes (Fig. 2A) revealed that this gene is absent in four species: Ciona intestinalis, Ciona savignyi, Drosophila melanogaster and Saccharomyces cerevisiae. Interestingly, a similar analysis conducted by us (Figs. S1–S3) for some other ALS risk factors such as angiogenin (ANG), Interferon Kappa (IFNκ) and NADH:Ubiquinone Oxidoreductase Subunit B6; (NDUFB6) reveals that these same genes are also absent in the four species mentioned above. At present it is unclear as to the significance of these findings.

It is worthwhile to note that C9orf72 does not have any paralogs in vertebrates and there is no evidence of gene family expansion or contraction events that have affected this gene. The ρ-value for the C9orf72 gene family was calculated to be 0.9410, suggesting that this gene has not gone through any significant gene gain or gene loss event during the course of vertebrate evolution. C9orf72 appears to be an ancient gene and the presence of an ortholog in Caenorhabditis elegans suggests that it must have been present in the last common ancestor of all vertebrates. It also appears that C9orf72 originated early in eukaryotic evolution since most extant eukaryotes possess only a single copy of the gene.

Statistically significant similarity between two protein sequences is frequently used to infer homology and therefore a common functional role. From the PSI-Blast (Altschul et al., 1997) results, all primate C9orf72 protein sequences have a pairwise identity of 99% with the human sequence and the mammals listed have alignment scores between 98% and 96% (Fig. 2B). The marsupials listed have scores of 94–96%, the birds have 93% and reptile proteins score 92%. Xenopus laevis scores 85% and the fish have scores ranging between 72–77%. This relatively high degree of protein similarity between species that are evolutionarily distant suggests there is a high degree of functional constraint on the protein. The residues are conserved all through the length of the protein suggesting that C9orf2 functions as a single domain protein. C9orf72 isoform containing 481 amino acids seems to be the major isoform expressed across all species. Given the high degree of sequence homology at the protein level and the consequent functional restraint on the protein, we wanted to know if this homology is reflected at the genetic level as well. We started by looking for conservation at the chromosome level followed by comparison of the different elements of the gene.

Synteny

Chromosome level comparison of the genes/markers between human and mouse revealed that almost 50% of the genes/markers were common to both species. Cinteny (Sinha & Meller, 2007) identified 19 syntenic blocks and computed a reversal distance of 5 (the minimum number of reversals to translate from one genome to another). Analysis of synteny around the C9orf72 loci showed that the gene is present in a conserved region of size 1.4 Mb in human (chromosome 9), whereas its mouse orthologue is present in a conserved region of size 2.1 Mb (chromosome 4).

The Synteny Database (Catchen, Conery & Postlethwait, 2009) was also used to identify syntenic regions between human and mouse (Fig. 3A). The synteny trace built by the software covered a region of 7.4 Mb in human (chromosome 9) and 7.9 Mb in mouse (chromosome 4). Concomitant with the larger syntenic region, the software detected an orthologous pairwise cluster comprising of 39 gene pairs (Table S1) that are conserved between the source (human) and the outgroup (mouse) genomes. The gene homology matrix (Fig. 3B) showed collinear regions (diagonal lines yellow squares) representing homologous genes between the human and mouse chromosomes. The gaps in diagonal regions (shown as blue squares) correspond to insertions that have arisen either through translocation or deletions of genes. Neither the two programs used above nor Ensembl has an annotated genome for Fugu. In order to compensate for the lack of an annotated Fugu genome, we used Ensembl to calculate syntenic regions from pairwise whole genome alignments of the human genome with stickleback, Tetraodon and zebrafish genomes. The syntenic regions calculated using Ensembl are in agreement with those detected by Cinteny and The Synteny Database.

Alternative splicing

The human C9orf72, through alternative splicing is transcribes five mRNA splice variants (Fig. 1). Of these, three transcripts are protein-coding variants with a validated polyadenylation tail at the 3′ end. The presence of the long and short isoforms of C9orf72 has been validated using isoform-specific antibodies (Xiao et al., 2015; Davidson et al., 2017). Improved detection using these antibodies helped reveal distinct subcellular localisation of the two protein isoforms as well as showed that they had distinct biochemical profiles (Xiao et al., 2015; Davidson et al., 2017). The remaining two spliced variants are processed transcripts but without an annotated ORF. It is currently assumed that these transcripts do not translate into viable protein products (DeJesus-Hernandez et al., 2011; Renton et al., 2011).

The existence of different C9orf72 splice variants of course will need to be experimentally validated by testing these antibodies on C9orf72 knock-in or knock-out tissues. In mouse, there are seven mRNA splice variants of which three are predicted to be protein coding and four not translated. Amongst these four variants of the mouse, two are predicted to undergo nonsense-mediated decay, as they are believed to retain intronic sequence relative to other coding variants. Alternative splicing of the Fugu gene is predicted to produce four transcription variants all of which are presumed to be translated to protein products (Fig. 1).

The longest mouse and human transcripts have a pairwise identity of 78%, the human vs. Fugu scored 70% and mouse vs. Fugu scored 74%. In human transcript 1, the start site is before exon 1 whereas in transcripts 2 and 3 it is found following exon 1—resulting in alternatively spliced transcripts that vary in their sequences (Fig. 1). The hexanucleotide repeat, GGGGCC, the expansion of which has been shown to be the pathologic link between ALS and FTD, is located in the core promoter region of human C9orf72 splice variant 1 and in intron 1 of transcript variants 2 and 3. In mouse, all of the variation is found prior to exon 3 apart from one, short, additional sequence in the final exon of transcript 3. In Fugu, there is little sequence and structural variation. Transcripts 1 and 2 are almost identical to each other except for the addition of a short sequence between exons 1 and 2 in transcript 1. Exon 1 in transcripts 2 and 4 are quite different in Fugu. Other than that, there are no differences in these two transcripts. Transcript 3 is the most unique in Fugu (Fig. 1).

Conservation of intron 1

The disease causing hexanucleotide repeat expansion occurs in intron one of C9orf72 in humans. Since it has been suggested that the repeat expansion may affect gene expression of C9orf72 (DeJesus-Hernandez et al., 2011; Renton et al., 2011) it is necessary to establish whether there is any conservation of intron 1 sequences as this would have implications when creating animal models.

We compared the intron 1 sequences from different species. High levels of sequence conservation in noncoding DNA sequences compared between species can be interpreted as evidence for functional constraints. We found that there are no clear conserved sequences generated by aligning the intron sequences of all the species to each other simultaneously. Because of this we started initially by looking at basewise conservation amongst the primates to evaluate potential signatures of evolutionary selection at particular positions of nucleotide bases. The alignment showed that significant conservation (blue bar graph) has been maintained throughout the sequence, as can be seen in Fig. 4. Also seen are sites that are fast evolving (red graph) but these are by far quite few in number and are present in mainly the intronic regions of the sequence. When this alignment was extended to the vertebrates, we found that the conserved residues were now confined mainly to the exons.

When comparing vertebrates, the introns seem to have undergone faster evolution than would be expected under neutral drift (evolutionary change that is caused not by natural selection but due to genetic drift of neutral mutant alleles). A region of around 130 bp in intron 1 appears to have higher conservation than the rest of the sequence. A closer look at this region showed a transcription factor match to ETS1. The Fugu sequence was then aligned with just the human sequence to try and find key, ancestral, conserved regions, if any. Pairwise global alignment revealed that less than 10% of the bases in the shorter Fugu sequence align with longer human intron 1 sequence (Fig. S4). The alignment improved when the sequences were analysed using local alignment tools, looking for local similarities instead of global alignments (Fig. S5). A similar pattern was observed when the human and mouse intron 1 sequences were aligned (Figs. S6 and S7).

Next we looked at the region surrounding the hexanucleotide repeat unit. Human intron 1 sequence was truncated to include only the region surrounding the repeat unit sequence. 200 bases both upstream and downstream of the repeat unit, plus the hexanucleotide sequence itself were used rather than the entire 6,622 bp sequence. The repeat region was then locally aligned to the entire Fugu intron sequence (576 bp). This alignment showed that just before the start of the repeat sequence (from position 201) there is a region of 24 nucleotides that has around 67% sequence identity (Fig. 5A). Repetitive motifs/elements are normally masked prior to sequence comparison to prevent spurious alignments. Since the Fugu sequence does not contain the repeat expansions, the hexanucleotide repeat units were removed from the human sequence and the two sequences were re-aligned. The removal of the repeat unit resulted in a sequence of 40 nucleotides on either side of the original repeat site with about 63% sequence identity (Fig. 5B). When the shortened version of the human sequence and the mouse intron sequences were locally aligned, we identified a 68 bp region with 67% sequence identity, which spans the repeat site (Fig. 5C).

The hexanucleotide repeat expansions in the non-coding region of C9orf72 gene are believed to be causal for ALS and FTD. So, we compared the presence/absence of this signature motif in the C9orf72 gene across species. The C9orf72 hexanucleotide GGGGCC repeat unit is found 241 bp into the sequence of intron 1 in human transcript variants 2 and 3, but not transcript 1. The intron sequence from transcript 2 is the longest, and was therefore used in our analyses. An unrepeated GGGGCC was also found further downstream in the sequence of human transcript 2 at position 6,245. The GGGGCC sequence is also found in the intron 1 sequence for chimpanzee, human’s closest relative, at position 2,220. However, the repeat unit is only found once and this position is significantly more downstream compared to the human sequence. Surprisingly, the promoter region (before exon 1) in chimpanzee bears three of the hexanucleotide repeat units occurring in tandem. The next closest species to human, gorilla, does not contain the hexanucleotide sequence in intron 1 although it is found in other primates including macaques, orangutans and marmosets. Interestingly, mouse, with the longest intron 1 amongst all the organisms compared, has two GGGGCC repeat units in this region, albeit not in tandem. The GGGGCC sequence was not seen in intron 1 of any of the Fugu transcripts.

Identification of putative immediate upstream elements

In the absence of biological data from experimental methods such as protein-binding microarrays, Chromatin Immunoprecipitation followed by Deep sequencing (ChIp-seq), we took a computational approach to identify the binding sites of DNA-associated proteins using DiAlign TF and MatInspector (Cartharius et al., 2005). To this end, the human C9orf72 gene up to the end of intron 1 together with 2,000 bp of the 5′-flanking region was aligned with the corresponding murine and Fugu sequences. The lack of sequence conservation throughout the alignment resulted in only one TF-binding site (within the aligned regions) that was common to all the three sequences. The transcription factor found common in all the three sequences was ETS1. So, in an attempt to find more transcription factor binding sites that are common to all three sequences, the sequences were scanned by MatInspector without performing any alignment. The matrices used for analysis were filtered based on the tissues they are associated with. The tissues chosen were brain, nervous system, neuroglia, neurons and spinal cord. A core similarity of 1.0 was used for the search. The program found a total of 46 TFs common to all three sequences out of which 18 general regulatory elements were omitted. Only those likely to be preferentially associated with nervous tissue-specific expression were retained (Table 3). The TFs identified play an important role in the development of the nervous system including neuronal differentiation, migration and fate specification (Table 3).

We also used MEME (Bailey & Elkan, 1994), Tomtom (Gupta et al., 2007) and GOMo (Buske et al., 2010), three motif-based sequence analysis tools from the MEME suite of programs to analyse the upstream regulatory region in human, mouse and Fugu. MEME (Bailey & Elkan, 1994) identifies motifs based on the assumption that TF motifs are more likely to be conserved in the regulatory region of a set of orthologous genes. We identified five potential TF binding motifs located at similar positions within the alignment of all the three species (Fig. 6A). The motifs discovered by MEME (Bailey & Elkan, 1994) (Fig. 6B) were searched against databases of known TF motifs using Tomtom (Gupta et al., 2007). TF matrix families identified by Tomtom (Gupta et al., 2007) (Table S2) are consistent with those determined by MatInspector (Cartharius et al., 2005). The MEME motifs were further analysed by the program GOMo (Buske et al., 2010) to determine significant/specific association of these motifs with genes that are linked to one or more gene ontology terms. This approach did not require the use of multiple sequence alignments as each comparative sequence was independently queried for association between a TF and a gene ontology (GO) term, thus avoiding drawbacks such as imperfect alignments or motif drift (Reddy, DeLisi & Shakhnovich, 2005). The GOMo analysis predicted several interesting biological roles for the motifs discovered by MEME. The significant predictions were neuron differentiation, ATP binding, axon guidance and negative regulation of transcription from RNA polymerase II promoter.

CpG islands and DNA methylation

CpG islands are typically common near transcription start sites and may be associated with promoter regions. These islands are rare in vertebrates because over time the methylated cytosines tend to spontaneously deaminate into thymines. Recent studies have shown that the CpG island, located at the 5′ end of the repeat expansion, is hypermethylated and could potentially lead to reduction of modulation of the disease phenotype (Xi et al., 2013; Belzil et al., 2013). It is believed that hypermethylation is neuroprotective because it inhibits transcription of the mutant C9orf72 leading to reduced RNA foci and dipeptide aggregate formation (Liu et al., 2014; McMillan et al., 2015; Xi et al., 2015). To this end, we looked at the CpG islands predicted by the UCSC genome browser (Human genome (version hg 38), Mouse genome (version mm 9), Tetraodon genome (version tetNig2), and Zebrafish genome (version danRer10)). What was instantly obvious was that only human C9orf72 gene has a CpG island, located at both 5′ and 3′ end of the region containing the repeat unit (Fig. 4). The 5′ CpG island seen in human intron 1 is absent in mouse, lamprey and zebrafish. To see if there was conservation of the CpG islands in primates in the promoter region, we did a pairwise alignment of the CpG island rich region at the 5′-end of human gene with that of that of the chimpanzee gene. We found that almost all the CpG dinucleotides that were investigated by Xi et al. (2013) are conserved in chimpanzee.