DNA RESEARCH 18, 177–187, (2011)
Advance Access Publication: 17 May 2011

doi:10.1093/dnares/dsr010

Exploration for Functional Nucleotide Sequence Candidates within
Coding Regions of Mammalian GenesRUMIKO Suzuki 1,2,† and NARUYA Saitou 1,2,*
Department of Genetics, School of Life Science, Graduate University for Advanced Studies, Mishima 411-8540,
Japan1 and Division of Population Genetics, National Institute of Genetics, Mishima 411-8540, Japan2*To whom correspondence should be addressed. Tel. þ81 559-81-6790. Fax. þ81 559-81-6789.
E-mail: saitounr@lab.nig.ac.jp
Edited by Katsumi Isono(Received 12 September 2010; accepted 14 April 2011)

1.

Introduction

The neutral theory of molecular evolution1,2 predicts that synonymous sites of codons are evolving
faster than non-synonymous sites because of thesmaller selective pressure. This is true in general;
however, several factors are known to influence on acertain region of a coding sequence and suppress
synonymous substitution.
One of the well-known factors is the codon bias
towards optimum codons. Optimum codons reflect
the composition of genomic tRNA pool.3 – 5 Because†

Present address: Department of Environmental and Preventive
Medicine, Oita University, Yufu 879-5593, Japan.

optimum codons are advantageous for fast and accurate translation, highly expressed or biologicallyimportant genes would have more optimum codons
than others.6 – 8 Changes from an optimum codon
to a non-optimum codon will be suppressed in suchgenes. Because optimum codons are similar inclosely related species, highly expressed genes tend
to show similar codon usage; therefore synonymous
substitution is lowered. In fact, the requirement fortranslational efficiency or accuracy enhances the
optimum codon usage and suppresses nucleotidechanges through purifying selection.5,6,9 – 11 Codon
bias towards optimum codons is strong in fastgrowing organisms such as Escherichia coli or
Saccharomyces cerevisiae, but generally weak in# The Author 2011. Published by Oxford University Press on behalf of Kazusa DNA Research Institute.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.
org/licenses/by-nc/2.5), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work
is properly cited.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

Abstract
The primary role of a protein coding gene is to encode amino acids. Therefore, synonymous sites ofcodons, which do not change the encoded amino acid, are regarded as evolving neutrally. However, if acertain region of a protein coding gene contains a functional nucleotide element (e.g. splicing signals),synonymous sites in the region may have selective pressure. The existence of such elements would be
detected by searching regions of low nucleotide substitution. We explored invariant nucleotide sequencesin 10 790 orthologous genes of six mammalian species (Homo sapiens, Macaca mulatta, Mus musculus,
Rattus norvegicus, Bos taurus, and Canis familiaris), and extracted 4150 sequences whose conservation
is significantly stronger than other regions of the gene and named them significantly conserved coding
sequences (SCCSs). SCCSs are observed in 2273 genes. The genes are mainly involved with development,
transcriptional regulation, and the neurons, and are expressed in the nervous system and the head andneck organs. No strong influence of conventional factors that affect synonymous substitution wasobserved in SCCSs. These results imply that SCCSs may have double function as nucleotide element
and protein coding sequence and retained in the course of mammalian evolution.
Key words: mammal; protein coding; nucleotide conservation

178

Exploration for Functional Nucleotide Sequence Candidates

2.

Materials and methods

2.1. Genome data
We obtained peptide and nucleotide sequences oforthologous genes of six mammalian species (Homosapiens, Macaca mulatta, Mus musculus, Rattus norvegicus, Bos taurus, and Canis familiaris) from Ensembl
database31 version 54 (http://May2009.archive.
ensembl.org/index.html). These species are selectedconsidering genome data quality and evolutionary
diversity. We selected one-to-one type single copy
orthologs and compiled 10 790 orthologous genesets (Supplementary file S1). We then constructed
multiple alignments of peptide sequences using
ClustalW32 and constructed nucleotide alignments

based on the peptide alignments. From the nucleotide alignments, we extracted sequences that areinvariant among all the species.
2.2. Identification of SCCSs
We performed permutation simulation to identify
SCCSs, or abbreviated as SCCSs, which are invariant
longer than 10 codons. This length is set to confine
the permutation run time within a feasible range.
For an N-codon long alignment, we generated a
non-redundant series of random numbers from 1 toN, and permuted codon columns (rows of codons inthe same site of the alignment) according to the generated random numbers. In this process, gap sites are
fixed and the rest of the sites are permuted. Then thelength and numbers of invariant sequences in the
permuted alignment are counted. We repeated thisprocess 500 000 times per ortholog set and took
averages of the frequency of invariant sequences. We
used the length and averaged frequency of invariant
sequences obtained from the permutation results as
random expectation, and evaluated the probabilityof invariant sequences in the original alignment
based on the expectation. This approach helps identify sequences whose conservation is rare to occur inthe substitution background of each alignment.
Multiple testing correction of p-values is done by
FDR (false discovery rate).33 Then we identified invariant sequences longer than 10 codons and P , 0.01
as SCCSs.
2.3. Analysis on GO, InterPro, and KEGG pathways
Protein coding genes that contain at least one SCCS
is named SCCS genes and those that do not containan SCCS is named non-SCCS genes. We used Fatigo
webservice(http://babelomics.bioinfo.cipf.es/functional.html) to identify GO terms, InterPro
codes, and KEGG pathways that are significantly
enriched with the SCCS gene group or with the nonSCCS gene group. Fatigo accepts a list of Ensembl
gene IDs as input and provides P-values for enrichment of the above terms. P-values are calculated by
Fisher’s exact test and corrected by FDR. We used
Ensembl gene IDs of H. sapiens as input and performed two-tailed comparison between SCCS genes
and non-SCCS genes.
2.4.

Analysis on preferred codons and average codon
degeneracy of SCCSs
We defined preferred codons as the most frequently
used codons for a given amino acid referring to Codon
Usage Database (http://www.kazusa.or.jp/codon/
index.html) provided by Kazusa DNA ResearchInstitute. Because the codon usage pattern was
similar among the six species we used, the codon

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

species with slow growth rate small population
size.12,13
Another factor is exonic splicing enhancer or silencer, which are splicing signals embedded inexons.14,15 Existence of such elements lowers thesynonymous substitution.16,17 In addition, ultraconserved elements (UCEs), which majorly reside in
non-protein coding regions, sometimes extend tocoding regions.18 In mammals, UCEs are reported toexist near to or overlap with genes associated with
nucleotide binding, transcriptional regulation, RNArecognition motif, zinc finger domain, and homeobox
domain.18 – 20 Hox genes also contain long conserved
nucleotide regions other than UCEs outside thehomeobox domain.21
Although the primary role of protein coding regionis to encode amino acids, there may be also functional
nucleotide elements embedded within coding
regions. For example, transcription-factor-bindingsites are found in coding regions,22 messenger RNAs
are targeted by various post-transcriptional regulations,23 and the requirement for a specific secondary structure for RNA editing decreases synonymous
substitution.24,25
Functional nucleotide elements are extensively
explored in the non-coding regions,26 – 30 but lessstudies have been done to explore probable functional
elements within the coding regions. We extracted significantly conserved coding sequences (SCCSs) fromorthologous genes of six mammalian species(human, rhesus macaque, mouse, rat, cow, and dog),
and compared genes containing SCCSs and genes
without SCCSs. Analyses on gene ontology (GO),InterPro codes, and KEGG pathways enlighten difference between the two gene groups. We also investigated RNA secondary structures, codon preference,
GC content, exonic splicing signals, and gene
expression of SCCSs to survey the influence of thesefactors.

[Vol. 18

No. 3]

R. Suzuki and N. Saitou

table of human was used as the representative. We
counted the number of preferred codons in a SCCS
and divided it by the codon length of the SCCS, andthen used the quotient as the ratio of preferredcodons. The average codon degeneracy is calculated
by summing up the degeneracy of each codon anddividing it by the codon length of the SCCS.

2.6. Evaluation of exonic splicing enhancer density
We obtained 238 hexamers from RESCUE-ESE Web
Server34 as candidates of exonic splicing enhancers.
We counted the number of the hexamers in the
SCCSs and the rest of the coding regions of the
10 790 genes (human sequences). We also measured
the total nucleotide length of the SCCS and the otherregions and applied the chi-square test for the ratio of
hexamers.

to count the length and frequency of invariant
sequences after random change of loci and use theresult as the random expectation.
We used the frequency distribution of invariant
sequences in the permuted alignments as therandom expectation and evaluated probability ofinvariant sequences in the original alignments, anddefined invariant sequences longer than 10 codons
and whose probability is below 0.01 (corrected by
FDR) as SCCSs.
In total, 4150 SCCSs (192 306 bp) were obtained
from 2273 alignments of 10 790 orthologous genesets (Supplementary file S2). This occupies 0.94% of
the coding region of the 10 790 genes. Table 1
shows the number of SCCSs per gene and thenumber of genes that contain that number of SCCSs.
Figure 1 is a graph of lengths and numbers of SCCSs(grey bars). Black dots indicate the random expectation obtained from the permuted alignments.
In the permuted alignments, there are 141
sequences (their total length is 5550 bp) whose probability is below 0.01, which means the region size of
SCCS is 35-fold larger than this expectation(Supplementary file S3). The x2 test between SCCSs

Table 1. Number of SCCSs in coding genes

SCCS genes

2.7. Analysis on gene expression
We used EGenetics (http://www.nhmrc.gov.au/
your_health/egenetics/index.htm) to investigate
gene expression of SCCS genes and non-SCCS genes.
Human anatomical system data, which give information about in which organs a gene is expressed,
were obtained from EGenetics database by Ensemble
Biomart. We counted how many of SCCS genes andnon-SCCS genes are expressed in each organ anddivided the numbers by the total number of SCCS
genes and non-SCCS genes, respectively. Then we performed the Fisher’s exact test to evaluate whether thedifference between SCCS and non-SCCS genes is significant. All P-values were corrected by FDR.
3.

Non-SCCS genes

Number of SCCSs(per gene)

Number of genes

1
2
3
4
5
.5

1366
475
219
105
42
66
8517

Results

3.1. Identification of SCCSs
If an alignment has a high ratio of conservation,
long invariant sequences may occur easily, and viceversa. Therefore, the rareness of invariant sequences
differs depending on the conservation backgroundin each alignment. To compensate this, we performed
permutation simulation. The idea of permutation is

Figure 1. The length and number of SCCSs. X and Y-axes represent
the length and frequency of SCCSs, respectively. Grey bars show
the length and the number of SCCSs. Black dots indicate thenumber of sequences with probability below 0.01 in the
permuted alignments.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

2.5. Prediction of RNA secondary structures
We computationally predicted secondary structuresand free folding energy of SCCSs using Vienna RNA
software package34 (http://www.tbi.univie.ac.at/~ivo/RNA/). Because folding free energy varies
depending on the sequence length, we constructed
free energy distribution by 1000 randomly chosensequences for each length (33 – 246 nucleotides).
The P-value for a given free energy was evaluatedbased on these distributions. Multi testing correction
is done by FDR.

179

180

Exploration for Functional Nucleotide Sequence Candidates

and the permutation result showed a significant
difference (P , 2.2E216).
If invariant and variant sites distributed randomly inthe original alignment, the frequency of invariant
sequences in the permuted alignments would
show similar frequency to the original alignment
because of the randomness at the start point. The
difference before and after the permutationsuggests that the distribution of invariant sites in theoriginal alignments is rather clustered than being
random.

[Vol. 18

3.2.

GO, InterPro, and KEGG pathways enriched in
SCCS containing genes
We used Fatigo web service to investigate differencein GO, InterPro codes, and KEGG pathways between
the SCCS genes and non-SCCS genes. The difference
was evaluated by the two-tailed Fisher’s test, and Pvalues were corrected by FDR. Tables 2 and 3 show
GO terms, InterPro codes, and KEGG pathways significantly enriched (P , 0.01) in SCCS genes.
The all terms in GO Biological process section isrelated to developmental process. One term of

Table 2. GO terms significantly (P , 0.01) enriched with SCCS genes
Terms

P*

In SCCS containing
genes (%)aIn non-SCCS containing
genes (%)b

Foldc

Biological process
7.75E205

0.53

0.04

13.25

2.50E204

0.48

0.04

12.00

GO:0060070: Canonical Wnt receptor signalling pathway

2.11E203

0.4

0.04

10.00

GO:0035137: Hindlimb morphogenesis

4.88E204

0.53

0.06

8.83

GO:0001702: Gastrulation with mouth forming second

1.78E203

0.44

0.05

8.80

GO:0009954: Proximal/distal pattern formation

3.80E204

0.57

0.07

8.14

GO:0031128: Developmental induction

1.03E203

0.53

0.07

7.57

GO:0048593: Camera-type eye morphogenesis

1.22E205

0.88

0.13

6.77

GO:0021510: Spinal cord development

2.37E204

0.75

0.13

5.77

GO:0031016: Pancreas development

1.85E203

0.62

0.12

5.17

GO:0014704: Intercalated disc8.43E203

0.35

0.04

8.75

GO:0043198: Dendritic shaft

2.38E203

0.48

0.06

8.00

GO:0030425: Dendrite

2.38E203

2.29

1.1

2.08

GO:0043025: Neuronal cell body

2.38E203

2.24

1.11

2.02

GO:0015629: Actin cytoskeleton

1.15E203

3.39

1.8

1.88

GO:0043005: Neuron projection

1.15E203

4.09

2.32

1.76

GO:0035254: Glutamate receptor binding

3.41E203

0.35

0.02

17.50

GO:0005072: Transforming growth factor beta receptor,
cytoplasmic mediator activity

8.92E203

0.31

0.02

15.50

GO:0004843: Ubiquitin-specific protease activity

3.71E203

0.48

0.07

6.86

GO:0031625: Ubiquitin protein ligase binding

6.44E203

0.62

0.14

4.43

GO:0003725: Double-stranded RNA binding

6.44E203

0.62

0.14

4.43

GO:0042054: Histone methyltransferase activity

8.64E203

0.57

0.13

4.38

GO:0050825: Ice binding

6.50E211

2.86

0.76

3.76

GO:0004221: Ubiquitin thiolesterase activity

4.18E204

1.01

0.27

3.74

GO:0005199: Structural constituent of cell wall

1.34E203

0.92

0.25

3.68

GO:0003682: Chromatin binding

4.08E208

2.15

0.59

3.64

Cellular components

Molecular functiona

The percentages of SCCS genes that have the GO term.
The percentages of non-SCCS genes that have the GO term.
c
The fold of ‘a’ to ‘b’. Terms are listed in the descending order of the fold difference. Terms with the highest 10-folds areshown for Biological Process and Molecular function.
*Probability for enrichment of the GO term in the SCCS group.
b

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

GO:0035136: Forelimb morphogenesis
GO:0035115: Embryonic forelimb morphogenesis

No. 3]

R. Suzuki and N. Saitou

181

Table 3. InterPro and KEGG terms significantly (P , 0.01) enriched with SCCS genes
Terms

P*

In SCCS containing genes(%)a

In non-SCCS containing
genes (%)b

Foldc

InterPro
IPR003619: MAD homology 1, Dwarfin-type
IPR001827: Homeobox protein, antennapedia type,conserved site

1.55E204

0.44

0.01

44.00

1.77E203

0.44

0.04

11.00

IPR000569: HECT

1.26E205

0.75

0.07

10.71

IPR002077: Voltage-dependent calcium channel,
alpha-1 subunit

5.65E204

0.53

0.05

10.60

IPR010982: Lambda repressor-like, DNA-binding

1.49E203

0.48

0.05

9.60

IPR002343: Paraneoplastic encephalomyelitis antigen

9.91E205

0.7

0.08

8.75

IPR018359: Bromodomain, conserved site

9.75E204

0.57

0.07

8.14

IPR001487: Bromodomain

1.94E206

0.97

0.12

8.08

IPR017995: Homeobox protein, antennapedia type

7.36E204

0.62

0.08

7.75

IPR004088: K Homology, type 1

7.36E204

0.62

0.08

7.75

KEGG
1.05E203

0.92

0.24

3.83

6.75E203

0.88

0.29

3.03

hsa04520: Adherens junction

4.55E203

1.06

0.37

2.86

hsa05211: Renal cell carcinoma

9.88E203

0.92

0.33

2.79

hsa04120: Ubiquitin-mediated proteolysis

1.05E203

1.5

0.57

2.63

hsa04310: Wnt signalling pathway

2.22E204

2.07

0.79

2.62

hsa04360: Axon guidance

3.47E204

1.8

0.69

2.61

hsa04810: Regulation of actin cytoskeleton

1.07E203

2.07

0.94

2.20

hsa04010: MAPK signaling pathway

3.23E204

3.08

1.5

2.05

a

The percentages of SCCS genes that have the InterPro or KEGG term
The percentages of non-SCCS genes that have the InterPro or KEGG term.
c
The fold of ‘a’ to ‘b’. Terms are listed in the descending order of the fold difference.
*Probability for enrichment of the InterPro or KEGG terms in the SCCS group.
b

Cellular components (GO:0014704) mediate mechanical and electrochemical integration between cardiomyocytes and the rest of the five (GO:0043198,
GO:0030425, GO:0043025, GO:0015629, andGO:0043005) have an association with the neuron.
In the molecular function category, three terms(GO:0004843, GO:0031625, and GO:0004221)
are related to the ubiquitin system, two(GO:00042054 and GO:00003682) are associatedwith chromatin. Ubiquitins are known to be involved
not only with protein degeneration but also withsignal transduction, chromatin modification, and cell
cycle.
Of the ten Interpro codes listed in Table 3,
IPR001827 is related to ubiquitin and IPR002077
represents calcium channel and other eight are all
associated with DNA or RNA-binding functions thatmediate transcriptional regulation or chromatin
modification.
Tow KEGG pathways (hsa04340 and hsa04310) in
Table 3 are developmental signalling pathways,

hsa04120 is related to the ubiquitin system, andhsa04360 is involved in axon guidance, which wellcorresponds with the GO and Interpro terms.
Table 4 shows the terms that are significantly scarcein SCCS genes. In contrast to Tables 2 and 3, majorityof the terms are involved with metabolic processes.
Fatigo can also explore enrichment of micro RNAtarget and transcription-factor-binding sites; no significant item was found.
3.3. Overlap with UCEs
UCEs are defined as nucleotide sequences that are
absolutely conserved longer than 200 bp between
orthologous regions of the human, rat andmouse.18,19 UCEs are found in both coding andnon-coding regions. Precedent studies report thatgenes with low synonymous substitution or genes
overlapped with UCEs are associated with DNA
binding, RNA binding, transcription activity, andHomeobox.18,19 In our 10 790 genes, 4009 bp in

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011hsa03018: RNA degradation
hsa04340: Hedgehog signalling pathway

182

Exploration for Functional Nucleotide Sequence Candidates[Vol. 18

Table 4. GO, InterPro and KEGG terms significantly (P , 0.01) scarce in SCCS genes
Terms

P*

In SCCS containing
genes (%)aIn non-SCCS containing
genes (%)b

Foldc

Biological process
GO:0043039: tRNA aminoacylation

6.15E203

0.09

0.63

0.14

GO:0006725: Cellular aromatic compound metabolic process

7.31E203

0.48

1.3

0.37

GO:0051186: Cofactor metabolic process

6.05E203

0.66

1.59

0.42

GO:0055114: Oxidation-reduction process

5.86E204

2.11

3.87

0.55

GO:0019752: Carboxylic acid metabolic process

3.54E204

2.42

4.34

0.56

GO:0006082: Organic acid metabolic process

3.54E204

2.42

4.35

0.56

GO:0044255: Cellular lipid metabolic process

8.60E203

3.3

4.98

0.66

Cellular component
GO:0019866: Organelle inner membrane

2.38E203

0.66

1.72

0.38

GO:0044429: Mitochondrial part

1.27E203

1.63

3.24

0.50

GO:0005615: Extracellular space

2.87E203

2.29

3.93

0.58

GO:0001595: Angiotensin receptor activity

8.64E203

0

0.4

0.00

GO:0016616: Oxidoreductase activity, acting on the CH-OH
group of donors, NAD or NADP as acceptor

6.42E203

0.18

0.82

0.22

GO:0016614: Oxidoreductase activity, acting on CH-OH groupof donors

8.92E203

0.26

0.97

0.27

8.96E203

0

0.46

0.00

2.86E203

0.57

1.56

0.37

Molecular function

IPR002198: Short-chain dehydrogenase/reductase SDR
Kegg
hsa04060: Cytokine– cytokine receptor interactiona

The percentages of SCCS genes that have the GO, InterPro or KEGG term.
The percentages of non-SCCS genes that have the GO, InterPro or KEGG term.
c
The fold of ‘a’ to ‘b’. Terms are listed in the ascending order of the fold difference.
*Probability for enrichment of the GO, InterPro or KEGG terms in the SCCS group.
b

29 genes are found to overlap with UCEs. In the4009 bp region, 2835 bp in 22 genes overlap with
SCCSs (Supplementary file S4). Because SCCSs areconserved in six mammalian species, including thethree species referred for UCEs, the 2835 bp reflect
conservation in other three species (macaque, cow,and dog). We surveyed nucleotide sites in the
10 790 genes that are conserved among human, rat,
and mouse, and evaluated how many of them arealso conserved in macaque, cow, and dog. The
resulted ratio is 0.751; therefore the expected conservation is 3035 bp. This matches well with the observation. The regions overlapped with UCEs make only
1.47% of the entire SCCS regions. Other SCCSs
convey shorter but deeper conservation than UCEs.
3.4. SCCSs that form stable RNA secondary structures
There are cases that a secondary structure of mRNA
conveys functions.24,25 We examined secondary structures and free energy of the SCCSs using Vienna RNA
package. We found three SCCSs whose folding
energy were significantly low (Table 5).

Polg encodes a catalytic subunit of mitochondrial
DNA polymerase POLG. The POLG protein is the only
polymerase known to be involved in replication of
mtDNA. Gal3st3 encodes a member of the galactose-3-O-sulfotransferase protein family. This protein
exists on the membrane of Golgi apparatus. Smarcd3
encodes a protein of SWI/SNF family, whose
members display helicase and ATPase activities. This
protein is thought to regulate transcription of target
genes by altering the chromatin structure aroundthose genes.
Table 5. Genes containing SCCS with significantly (P , 0.001) low
free folding energy
Gene

LengthFreeenergy219.9

polg

DNA polymerase subunit gamma-1

36

gal3st3

Galactose-3-O-sulfotransferase 3

36

222.6

smarcd3

SWI/SNF-related matrix-associatedactin-dependent regulator ofchromatin subfamily D member

39

223.9

The gene names are represented by those of human.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

Interpro

No. 3]

R. Suzuki and N. Saitou

183

Table 6. Splicing signals in SCCS and non-SCCS regions of 10 790
genes

SCCS

Region size (bp)
192 306#Splicing signals
20 420Per nucleotide
0.106

Non-SCCS

20 566 520

2 183 544

0.106

Splicing signals in non-SCCS regions are counted on humansequences.

Figure 2 shows the probability density of the folding
free energy constructed by randomly extracted
sequences. Each line shows the free energy of thesequences of the same length as the above three
SCCSs. Gene names on the lines represent freeenergy of the SCCSs. This figure suggests that thesethree SCCSs have extremely low free energy and thatthese regions will form stable secondary structures.
These SCCSs may be conserved because of therequirement for the secondary structures. However,substitution restriction of this type might be modest
because a secondary structure can be retained byanother combination of nucleotides as far as the complementarity is maintained.

3.5.

The density of exonic splicing enhancers in SCCSs
and in the other coding regions
One of the well-known functional nucleotide
elements in the coding region is splicing signals. We
obtained 238 hexamers from RESCUE-ESE Web
server as candidates of exonic splicing enhancers,
and counted the number of hexamers in SCCSs andnon-SCCS regions of the 10 790 genes (Table 6).
Splicing signals in non-SCCS regions are counted onhuman sequences. We observed 20 420 hits ofsignals in 192 306 bp of SCCS regions and2 183 544 hits in 205 666 520 of non-SCCS
regions. The number of the signals per nucleotide isboth 0.106 for SCCS and non-SCCS, and there wasno significant difference (P ¼ 0.99, x2 test).

3.7.

Preferred codons, GC content, and codon
degeneracy of SCCSs
Codon usage biases towards optimum codons are
known to suppress synonymous substitution.
Optimum codons reflect the composition of the
genomic tRNA pool and are advantageous for translation efficiency or accuracy. As the approximate
index of optimum codons, we used preferredcodons, or most frequently used codon for an aminoacid, and evaluated the preferred codon fraction in
SCCSs. In SCCS regions, 28 188 of 64 102 codons
are preferred codons and in non-SCCS regions,
2 961 482 of 6 855 505 codons are preferredcodons. The ratio of preferred codons in SCCS andnon-SCCS regions are 0.440 and 0.432, respectively.
The difference is significant at the 0.05 significancelevel (P ¼ 0.013) but the difference of the ratios ismerely 0.008. We also observed that the ratio of preferred codon decreases as the length of SCCS increases(Fig. 3). Judging from these results, SCCSs are unlikely
to be retained solely by codon preference.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

Figure 2. The probability density of the folding free energy. Thisgraph shows the probability density of the folding free energy
for 36 and 39 nucleotide-long sequences. The three SCCSs
with significantly low free energy are indicated on the graph.
Probability density was created using software package R.38

3.6. Gene expression
We investigated the difference of gene expression
between SCCS genes and non-SCCS genes referring
to anatomical system data of EGenetics, which givequalitative information about in what organs a geneis expressed. We counted the number of SCCS genes
and non-SCCS genes expressed in the organs and performed the Fisher’s exact test as described in Materials
and method section.
We compared the percentages of genes that are
expressed in each organ. Table 7 shows organs inwhich significantly higher percentage of SCCS genes
are expressed compared with non-SCCS genes. In
general, SCCS genes are expressed in a wider varietyof organs. This observation agrees with a previous
study.35 Only in medulla oblongata and trophoblast,
non-SCCS genes showed significantly higher percentage than SCCS genes (data not shown).
Seven organs (amygdala, spinal cord, cerebellum
cortex, cerebellum, frontal lobe, pituitary grand, andsympathetic chain) in Table 7 are related with thenervous system and six organs (cochlea, trabecularmeshwork, hypopharynx, larynx, tongue, andthyroid) are associated with head and neck.

184

Exploration for Functional Nucleotide Sequence Candidates[Vol. 18

Table 7. Organs in which significantly (P , 0.001) higher percentage of SCCS genes are expressed compared with non-SCCS genes
Organ

P*

SCCS

Foldc

Non-SCCSab#Expressed

Percentage#Expressed

Amygdala

1.32E211

201

9.86

318

Percentage
5.37

1.84

Cochlea

1.24E222

436

21.39

723

12.21

1.75

Small intestine

5.28E203

49

2.40

86

1.45

1.66

Amnion1.03E204

102

5.00

183

3.09

1.62

Amniotic fluid

9.12E207

175

8.59

321

5.42

1.58

Spinal cord

7.21E205

129

6.33

243

4.10

1.54

Artery9.80E205

133

6.53

254

4.29

1.52

Cerebellum cortex

2.74E203

83

4.07

160

2.70

1.51

Cerebellum

3.62E208

277

13.59

543

9.17

1.48

Trabecular meshwork

1.38E205

199

9.76

399

6.74

1.45

Frontal lobe

1.63E228

899

44.11

1803

30.45

1.45

Hypopharynx

1.75E206

246

12.07

497

8.39

1.44

1.67E209

421

20.66

877

14.81

1.39

9.65E206

271

13.30

575

9.71

1.37

Breast

5.95E230

1133

55.59

2430

41.04

1.35

Larynx9.05E213

642

31.50

1385

23.39

1.35

Tongue

4.55E207

377

18.50

816

13.78

1.34

Smooth muscle

1.36E205

307

15.06

670

11.32

1.33

Thyroid

3.36E223

1076

52.80

2375

40.11

1.32

Adrenal gland

4.51E206

362

17.76

801

13.53

1.31

a

The percentages of SCCS genes that are expressed in the organ.
The percentages of non-SCCS genes that are expressed in the organ.
c
The fold of ‘a’ to ‘b’ Terms are listed in the descending order of the fold difference.
*Probability for enrichment of the expressed genes in the SCCS group.
b

Figure 3. The ratio of preferred codons in SCCSs. X-axis represents
the length of SCCSs, and Y-axis represents the ratio ofpreferred codon of the sequences. Classes whose sample size
, 20 were combined. Error bars represent 1 SE.

In mammals, of the influence GC content on
nucleotide change as a result of CpG hyper mutability.
We investigated GC content of SCCSs. GC content in
the first (GC1), second (GC2), and third position(GC3) of codons show different patterns along thesequence length (Fig. 4). GC1 is mostly constant but
GC2 increases while GC3 decreases as the length of
SCCS increases. Because mammalian genomes prefer
GC-ending codons, the decrease of GC3 corresponds
to the decrease of preferred codons. The decrease of
GC3 seems to be complementary with the increaseof GC2 because GC content as a whole is constant(Supplementary Fig. S1).
Conservation of SCCSs may occur by chance in aregion where amino acid constraint is strong andcodon degeneracy is low. We investigated codon degeneracy of SCCSs to examine this possibility. The averagedegeneracy is between three and four and increases
as the sequence length increases (Fig. 5). This result
suggests that even if the first and second positions ofcodons are restricted, the third position has enoughfreedom to change. Therefore, it is unlikely that SCCSs
are conserved because of the amino acid constraint
combined with the low degeneracy of codons.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

Pituitary glandSympathetic chain

No. 3]

R. Suzuki and N. Saitou

185

Figure 4. GC content of the first (GC1), the second (GC2), and thethird (GC3) position of codons in SCCSs: (A) GC1, (B) GC2, (C)GC3. X-axis represents the length of SCCSs, and Y-axis represents
GC content of the sequences. Classes whose sample size , 20
were combined. Error bars represent 1 SE.

4.

Discussion

GO terms, InterPro codes, and KEGG pathways
enriched with SCCS genes show a strong commitment
to the developmental process, transcriptional regulation, and the neurons. Genes associated with transcriptional regulation or the neurons are known tohave a low synonymous substitution ratio. Thisphenomenon is discussed in relation with codon

biases to improve translational efficiency or accuracy.
However, our analysis on the preferred codons in
SCCSs suggests that codon preference is not likely
the major factor influencing on the conservation of
SCCSs.
Analyses on the ratio of preferred codons, GC
content, and codon degeneracy enlighten the characteristics of SCCSs. The ratio of preferred codons
decreases as the SCCS length increases. Drummondand Wilke36 investigated the correlation between
synonymous substitution rate (dS) and fraction of
optimum (Fop) codons, and detected negative correlations between dS and Fop in rodents (mouse andrat) and positive correlation in human – dog comparison. If the SCCSs have the same trend as the rodents ofthe previous study, the ratio of preferred codons in
SCCSs should be high; however, our result is not.
There was no factor that would lower nucleotide substitution in GC content and codon degeneracy.
Methodological difference is that our research
focused on local and complete conservation of
nucleotides instead of the dS in the entire region ofa gene and that we investigated conservation among
the six mammalian species instead of a pair-wisecomparison. The difference of results may suggestthat factors underlying local and strong conservationsuch as SCCS differ from the factors working on thegene-wide conservation.
Makalowski et al.37 showed a correlation between
synonymous substitution rate (dS) and non-synonymous substitution rate (dN). Such correlations may
occur when the constraint on a certain nucleotide
sequence is so strong that dN is also lowered.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

Figure 5. Codon degeneracy of SCCSs. X-axis represents thelength of SCCSs and Y-axis represents the averaged codon
degeneracy of the sequences. Classes whose sample size , 20
were combined. Error bars represent 1 SE.

186

Exploration for Functional Nucleotide Sequence Candidates

The usage of relatively rare codons and strong local
conservation of SCCSs may be preferable as regulatory
signals. The fraction of SCCSs in the coding region of
the 10 790 genes is 0.94%. This fraction is so smallthat it would not have an influence on conventional
evolutionary analysis. Although the fraction is small,or because the fraction is small, SCCSs may havepotential as regulatory elements.
Acknowledgements: We thank Drs Kenta
Sumiyama, Kiyoshi Ezawa, and Hiroshi Akashi for
their insightful suggestion and discussion.
Supplementary data: Supplementary Data are
available at www.dnaresearch.oxfordjournals.org.

Funding

References
1. Kimura, M. 1983, The Neutral Theory of Molecular
Evolution. Cambridge University Press: Cambridge.
2. Nei, M. 1987, Molecular Evolutionary Genetics. Columbia
Univ. Press: New York.
3. Ikemura, T. 1982, Correlation between the abundanceof yeast transfer RNAs and the occurrence of therespective codons in protein genes. Differences insynonymous codon choice patterns of yeast andEscherichia coli with reference to the abundance of isoaccepting transfer RNAs, J. Mol. Biol., 158, 573–97.
4. Ikemura, T. 1981, Correlation between the abundanceof Escherichia coli transfer RNAs and the occurrence ofthe respective codons in its protein genes: a proposal
for a synonymous codon choice that is optimal for the
E. coli translational system, J. Mol. Biol., 151, 389– 409.
5. Kanaya, S., Yamada, Y., Kinouchi, M., Kudo, Y., andIkemura, T. 2001, Codon usage and tRNA genes in eukaryotes: correlation of codon usage diversity with translation efficiency and with CG-dinucleotide usage asassessed by multivariate analysis, J. Mol. Evol., 53, 290–8.
6. Akashi, H. 1994, Synonymous codon usage in
Drosophila melanogaster: natural selection and translational accuracy, Genetics, 136, 927 –35.
7. Kurland, C.G. 1991, Codon bias and gene expression,
FEBS Lett., 285, 165 –9.
8. Hershberg, R., and Petrov, D.A. 2008, Selection on codon
bias, Annu. Rev. Genet., 42, 287–99.
9. Ikemura, T. 1985, Codon usage and tRNA content inunicellular and multicellular organisms, Mol. Biol. Evol.,2, 13– 34.

10. Sharp, P.M., and Li, W.H. 1987, The rate of synonymous
substitution in enterobacterial genes is inversely related
to codon usage bias, Mol. Biol. Evol., 4, 222– 30.
11. Akashi, H. 2003, Translational selection and yeast proteome evolution, Genetics, 164, 1291 –303.
12. Eyre-Walker, A.C. 1991, An analysis of codon usage inmammals: selection or mutation bias? J. Mol. Evol.,33, 442–9.
13. Sharp, P.M., Averof, M., Lloyd, A.T., Matassi, G., andPeden, J.F. 1995, DNA sequence evolution: the soundsof silence, Philos. Trans. R Soc. Lond. B Biol. Sci., 349,
241– 7.
14. Reed, R. 1996, Initial splice-site recognition and pairing
during pre-mRNA splicing, Curr. Opin. Genet. Dev., 6,
215– 20.
15. Blencowe, B.J. 2000, Exonic splicing enhancers: mechanism of action, diversity and role in human genetic diseases. Trends Biochem. Sci., 25, 106– 10.
16. Takahashi, A. 2009, Effect of exonic splicing regulation
on synonymous codon usage in alternatively spliced
exons of Dscam, BMC Evol. Biol., 9, 214.
17. Parmley, J.L., and Hurst, L.D. 2007, Exonic splicing regulatory elements skew synonymous codon usage nearintron-exon boundaries in mammals, Mol. Biol. Evol.,24, 1600– 3.
18. Bejerano, G., Pheasant, M., Makunin, I., et al. 2004,
Ultraconserved elements in the human genome,
Science, 304, 1321– 5.
19. Schattner, P., and Diekhans, M. 2006, Regions of
extreme synonymous codon selection in mammaliangenes, Nucleic Acids Res., 34, 1700– 10.
20. Lareau, L.F., Inada, M., Green, R.E., Wengrod, J.C., andBrenner, S.E. 2007, Unproductive splicing of SR genes
associated with highly conserved and ultraconserved
DNA elements, Nature, 446, 926 –9.
21. Lin, Z., Ma, H., and Nei, M. 2008, Ultraconserved coding
regions outside the homeobox of mammalian Hox
genes, BMC Evol. Biol., 8, 260.
22. Lampe, X., Samad, O.A., Guiguen, A., et al. 2008, An
ultraconserved Hox-Pbx responsive element resides in
the coding sequence of Hoxa2 and is active in rhombomere 4, Nucleic Acids Res., 36, 3214– 25.
23. Licatalosi, D.D., and Darnell, R.B. 2010, RNA processing
and its regulation: global insights into biological networks, Nat. Rev. Genet, 11, 75–87.
24. Bhalla, T., Rosenthal, J.J., Holmgren, M., and Reenan, R.
2004, Control of human potassium channel inactivation by editing of a small mRNA hairpin, Nat. Struct.
Mol. Biol., 11, 950 –6.
25. Delgado, M.D., Gutierrez, P., Richard, C., Cuadrado, M.A.,
Moreau-Gachelin, F., and Leon, J. 1998, Spi-1/PU.1
proto-oncogene induces opposite effects on monocyticand erythroid differentiation of K562 cells, Biochem.
Biophys. Res. Commun., 252, 383 –91.
26. Donehower, L.A., Slagle, B.L., Wilde, M., Darlington, G.,
and Butel, J.S. 1989, Identification of a conserved
sequence in the non-coding regions of many humangenes, Nucleic. Acids Res., 17, 699 –710.
27. Lipman, D.J. 1997, Making (anti)sense of non-coding
sequence conservation, Nucleic Acids Res., 25, 3580– 3.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011

This study was supported partly by Grant-in-Aid forscientific research from the Ministry of Education,
Culture, Sports, Science, and Technology of Japan toN.S.

[Vol. 18

No. 3]

R. Suzuki and N. Saitou

28. Levy, S., Hannenhalli, S., and Workman, C. 2001,
Enrichment of regulatory signals in conserved noncoding genomic sequence, Bioinformatics, 17,871 – 7.
29. Sandelin, A., Bailey, P., Bruce, S., et al. 2004, Arrays ofultraconserved non-coding regions span the loci of
key developmental genes in vertebrate genomes, BMC
Genomics, 5, 99.
30. Kryukov, G.V., Schmidt, S., and Sunyaev, S. 2005, Small
fitness effect of mutations in highly conserved noncoding regions, Hum. Mol. Genet., 14, 2221– 9.
31. Hubbard, T.J., Aken, B.L., Beal, K., et al. 2007, Ensembl
2007, Nucleic Acids Res., 35, D610 – 7.
32. Thompson, J.D., Higgins, D.G., and Gibson, T.J. 1994,
CLUSTAL W: improving the sensitivity of progressive
multiple sequence alignment through sequence
weighting, position-specific gap penalties and weightmatrix choice. Nucleic Acids Res., 22, 4673– 80.

187

33. Benjamini, Y., Drai, D., Elmer, G., Kafkafi, N., andGolani, I. 2001, Controlling the false discovery rate inbehavior genetics research, Behav. Brain Res., 125,
279– 84.
34. Fairbrother, W.G., Yeh, R.F., Sharp, P.A., and Burge, C.B.
2002, Predictive identification of exonic splicing enhancers in human genes, Science, 297, 1007 –13.
35. Plotkin, J.B., Robins, H., and Levine, A.J. 2004, Tissuespecific codon usage and the expression of humangenes, Proc. Natl. Acad. Sci. USA, 101, 12588– 91.
36. Drummond, D.A., and Wilke, C.O. 2008, Mistranslationinduced protein misfolding as a dominant constraint on
coding-sequence evolution, Cell, 134, 341– 52.
37. Makalowski, W., and Boguski, M.S. 1998, Synonymousand nonsynonymous substitution distances are correlated in mouse and rat genes, J. Mol. Evol., 47, 119 –21.
38. Ihaka, R., and Gentleman, R. 1996, R: a language for data
analysis and graphics, J. Comp. Graph. Stat., 5, 299–314.

Downloaded from dnaresearch.oxfordjournals.org at Idenkagaku Kenkyujo on August 9, 2011