Case Study 1: Construction of CNS Alignment, Including an SNP That Causes Human Disease

Based on coordinates of human SNP sites, dbCNS can construct multiple sequence alignments to evaluate evolutionary conservation of genomic regions, including specified sites. Aniridia (OMIM ID 106210) is a panocular disease characterized by a variable degree of iris/foveal hypoplasia, nystagmus, and ciliary body abnormalities. In a patient with aniridia and no exonic mutations or chromosomal abnormalities, direct sequencing of cis-regulatory elements active in various eye tissues revealed a single-nucleotide change in a conserved ocular enhancer, SIMO, located 150 kb downstream from PAX6 (Bhatia et al. 2013). The SNP that causes aniridia (Bhatia et al. 2013) is C>A at rs606231388 in dbSNP.

As we already showed in an example of “sequence extraction mode,” dbCNS extracted a 201-bp sequence, including this SNP site, from the reference human genome sequence (hg38) using “11:31664397>A” as a keyword (supplementary fig. S3B, Supplementary Material online). Using this sequence, output of the BLAST & alignment mode was generated with 38 gnathostome genomes (supplementary fig. S4A, Supplementary Material online). In this analysis, the “-num_alignments” option was set at two in order to count identified CNSs in each species. As a result, dbCNS identified at most one BLAST hit for each species (shown in [# of blast hits]) and automatically aligned them. The alignment showed that all BLAST hits of gnathostomes contain the PAX6 binding site and belong to the SIMO region (Bhatia et al. 2013), except for the partial sequence of Erpetoichthys calabaricus (reedfish). Then, we confirmed that these BLAST hits are identical to the human query sequence and form a CNS as a highly conserved part of the SIMO enhancer (Antosova et al. 2016).

The alignment (fig. 4A) confirmed that in this aniridia-related site, most tetrapods share the same nucleotide C and the mutation changed the human nucleotide from C>A. In addition, the alignment showed that all five snakes share A at this site. In this case, dbCNS can be used to detect CNS candidates with accelerated substitution rates. The estimated CNS tree (fig. 4B) suggested that in the snake lineage, branches leading to the common ancestor of the five snakes possessed an increased number of substitutions compared with peripheral branches. These findings imply that characteristics of the snake SIMO region were fixed before divergence of the major snake lineages. This CNS diversification in snake ancestors is consistent with their possible subterranean lifestyle (Da Silva et al. 2018) and the loss of opsins in the early stage of snake evolution (Simoes et al. 2015). In contrast, four subterranean mammals (species names are shown in red) showing convergent eye degeneration shared the nucleotide C with most other tetrapods (fig. 4A). In subterranean mammals, several CNSs near PAX6 loci and other transcription factors important for eye development exhibit accelerated substitution rates (Partha et al. 2017). In this analysis of the SIMO region, an accelerated substitution rate was suggested for the lineage leading to the subterranean mammal, Heterocephalus glaber (naked mole rat), compared with other eutherians (fig. 4B). The naked mole rat sequence is not placed next to related species probably due to its high sequence divergence. For more sophisticated analyses of accelerated substitution rates with user-defined tree topologies, users can employ state-of-the-art methods, such as RERconverge (Kowalczyk et al. 2019), using dbCNS outputs.

Fig. 4.

(A) Alignment of the main part of the SIMO region (supplementary fig. S4, Supplementary Material online). In the query sequence, YOURSEQ1, the SNP site is highlighted with a red background. (B) Phylogenetic tree based on sequences of the SIMO region (121 sites).

dbCNS produces a link to dbSNP (build 153) using BLAST-hit coordinates derived from the human genome (hg38). By clicking the link “11:31664297-31664497” located below “Human SNP in dbSNP:” in the output html file (supplementary fig. S4A, Supplementary Material online), it was confirmed that the aniridia-causing SNP site (rs606231388) is located in this human BLAST hit. Moreover, dbCNS can analyze SNPs identified in genome-wide association studies. For example, nasopharyngeal carcinoma-related SNP (Madelaine et al. 2018) can be analyzed using “3:169364845>A” (hg38) as a keyword (supplementary fig. S4B, Supplementary Material online).

Case Study 2: Detection of CNSs in Gnathostome Genomes

dbCNS can evaluate the existence or number of CNSs in genomes. In addition to the SIMO region, several CNSs were annotated as cis-regulatory elements that control expression of PAX6 in various tissues, including the eye. Bhatia et al. (2014) identified CNSs in the RCN1–PAX6 intergenic region by employing a strategy that analyzes gnathostome sequence conservation and tests identified CNSs of the elephant shark for enhancer activity using a combination of zebrafish and mouse transgenic studies.

Thus, we examined CNSs shared among other gnathostomes. Using 20 published CNS coordinates (supplementary table S1, Supplementary Material online), the dbCNS “Sequence extraction” mode reported CNSs from human genome data (supplementary table S2A, Supplementary Material online). The existence of identical CNSs was then evaluated for 38 gnathostome genomes using extracted CNSs as queries in the “BLAST & alignment” mode. The option “-num_alignments” was set at two to detect duplicated CNSs in each genome. A method for conducting multiple analyses and summarizing results is shown on the instruction page. BLAST hits were detected for all species analyzed only for agCNE2 and cre149, whereas all four teleost genomes lacked nine CNSs (fig. 5A). Although most BLAST hits were single, two hits were detected for several species, such as Podarcis muralis (common wall lizard), Equus caballus (horse), and Aotus nancymaae (Nancy Ma's night monkey).

Fig. 5.

Results of gnathostome analyses. (A) BLAST hits for CNS queries around human PAX6 locus. An arrowhead indicates the row of humans, sequences of which were used as queries. Phylogenetic positions of whole-genome duplications (VGD, vertebrate genome duplication; TGD, teleost genome duplication) follow Braasch and Postlethwait (2012). The heatmap was summarized by using a script available from the dbCNS instruction page. (B) Overview of CNS positions around PAX6 loci. The black line represents DNA. Red letters indicate CNS queries in humans and blue letters indicate CNS BLAST hits in nonhuman gnathostomes. Rectangles indicate the PAX6 locus (red) and adjacent RCN1 (green) and ELP4 loci (yellow). In P. muralis, thin horizontal arrows indicate putative duplicated regions. Arrows within gene loci indicate TSS.

These BLAST hits were mapped onto genomic regions of eight gnathostome species to determine the presence of CNSs around the PAX6 locus (fig. 5B). Summary statistics from those 20 analyses were generated by using our customized command-line scripts available from the dbCNS instruction page. In addition to their genomic positions, we identified CNSs by evaluating sequence alignments and bit scores. As a result, a single CNS was identified in this region, although P. muralis possessed six duplicated CNSs.

The six duplicated CNSs (agCNS9–13 and P2) of P. muralis formed a pair of blocks: an 11-kb region consisting of the six CNSs with the same order as in the human genome and a 37-kb region, including additional three CNSs (agCNS6–8) with reversed order. The arrangement of P. muralis CNSs appears to show vestiges of duplication and inversion, though CNSs of other gnathostomes were aligned in the same order as in the human genome. This conserved architecture shared among gnathostomes is probably important for the PAX6 system.

To illustrate the novelty of dbCNS, identified CNSs were compared with those estimated using a pioneering web tool in this field, mVISTA (Frazer et al. 2004; http://genome.lbl.gov/vista/mvista/submit.shtml), with special reference to the vestiges of duplication and inversion in the intergenic RCN1–PAX6 region. As far as we know, except for dbCNS, mVISTA is the only other web tool that can identify CNSs in response to user requests. Although mVISTA has a feature to identify novel CNSs, users must prepare sequences of interest for all species. When the intergenic RCN1–PAX6 region was compared among eight species used in gnathostome analyses (fig. 5B), mVISTA identified almost the same CNSs sets (supplementary fig. S6A, Supplementary Material online) as those identified by dbCNS (supplementary fig. S6B, Supplementary Material online). However, despite weak detection of two CNSs (agCNE7 in Lepisosteus and agCNE12 in Danio) not detected in dbCNS analyses, mVISTA could not identify 11 CNSs located within the 37-kb block of Podarcis. Vestiges of duplication and inversion prevented mVISTA from identifying these duplicated CNSs using multiple sequence alignments.

Case Study 3: Detection of Lineage-Specific CNSs from Teleost Genomes

dbCNS can detect lineage-specific CNSs. Due to the additional whole genome duplication in the teleost lineage (teleost genome duplication [TGD]) and its consequently increased rate of evolutionary divergence, teleost genomes lack many CNSs identifiable in other vertebrates (Lee et al. 2011). In fact, only three CNSs (agCNE9, agCNE13, and cre149 in fig. 5B) were identified around the PAX6b locus of the Danio rerio (zebrafish) genome in our gnathostome analysis. Using the keyword “PAX6b” for an analysis in the “Keyword search” mode, 164 CNSs conserved between zebrafish (D. rerio) and sticklebacks (Gasterosteus aculeatus) were listed. Among those, 30 zebrafish sequences (zs1–zs30 in supplementary table S2B, Supplementary Material online) had more than four BLAST hits when analyses were conducted in BLAST & alignment mode with our teleost taxon sampling (fig. 6A). Results of these 30 analyses are summarized on the right side of figure 6A.

Fig. 6.

Results of teleost analyses. (A) BLAST hits for CNS queries around the zebrafish PAX6b locus. All 30 zsCNSs were identified by comparing zebrafish and stickleback genomes in the ANCORA database (Engstrom et al. 2008). An arrowhead indicates the row of Danio rerio, sequences of which were used as queries. (B) Overview of CNS positions around zebrafish and medaka PAX6b and zebrafish PAX6a loci. Red letters indicate CNSs queries around zebrafish PAX6b loci and blue letters indicate CNS blast hits in other regions. (C) Ancestral PAX6 synteny blocks of teleosts. PAX6 synteny blocks were compared among bony vertebrate genomes using a conserved synteny browser, Genomicus ver. 98.01 (Muffato et al., 2010). Based on the PAX6 gene tree (Feiner et al. 2014), hypothetical ancestral states around the PAX6 locus were reconstructed using parsimony.

Single BLAST hits were detected in many cases (fig. 6A). When mapping BLAST hits of Oryzias latipes (medaka) chromosome 3 on the region around the PAX6b locus, 17 of 30 query CNSs of D. rerio had identical CNSs (blue letters in fig. 6B). However, two BLAST hits were detected when some of these 30 CNS queries were used, especially for D. rerio (fig. 6A). When mapping BLAST hits of D. rerio around the PAX6a locus in chromosome 25, 10 out of 30 query CNSs (blue letters in fig. 6B) had TGD-derived counterparts. The teleost PAX6a gene is known as the counterpart of the PAX6b gene derived from the TGD (Feiner et al. 2014). Given the preservation of the ancestral PAX6 synteny block around the teleost PAX6b locus (fig. 6C), and the singleton status of PAX6b-adjacent genes, RCN1 and ELP4, in the last common ancestor of teleosts (supplementary fig. S5, Supplementary Material online), counterparts of teleost RCN1 and ELP4 genes from TGD are considered lost from the region around the PAX6a locus (fig. 6C). This disappearance of adjacent genes, RCN1 and ELP4 counterparts, supports the hypothesis that in the D. rerio genome, these ten CNS counterparts function as regulatory elements of the PAX6a gene, as suggested by Kikuta et al. (2007).