Elucidation of phenotypic adaptations: Molecular
analyses of dim-light vision proteins in vertebrates
Shozo Yokoyama*†, Takashi Tada*, Huan Zhang‡, and Lyle Britt§

*Department of Biology, Emory University, Atlanta, GA 30322; ‡Department of Marine Sciences, University of Connecticut, Groton, CT 06340;
and §Alaska Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA 98195

Edited by Masatoshi Nei, Pennsylvania State University, University Park, PA, and approved July 14, 2008 (received for review March 12, 2008)

Vertebrate ancestors appeared in a uniform, shallow water envi-
ronment, but modern species fourish in highly variable niches. A
striking array of phenotypes exhibited by contemporary animals is
assumed to have evolved by accumulating a series of selectively
advantageous mutations. However, the experimental test of such
adaptive events at the molecular level is remarkably diffcult. One
testable phenotype, dim-light vision, is mediated by rhodopsins.
Here, we engineered 11 ancestral rhodopsins and show that those
in early ancestors absorbed light maximally (�max) at 500 nm, from
which contemporary rhodopsins with variable �maxs of 480 –525
nm evolved on at least 18 separate occasions. These highly envi-
ronment-specifc adaptations seem to have occurred largely by
amino acid replacements at 12 sites, and most of those at the
remaining 191 (�94%) sites have undergone neutral evolution.
The comparison between these results and those inferred by
commonly-used parsimony and Bayesian methods demonstrates
that statistical tests of positive selection can be misleading without
experimental support and that the molecular basis of spectral
tuning in rhodopsins should be elucidated by mutagenesis analy-
ses using ancestral pigments.

molecular adaptation � rhodopsin

The morphologies and lifestyles of animals in a wide range of environmental conditions have evolved to generate a striking
array of forms and patterns. It is generally assumed that these
variations have been driven by mutations, followed by positive
Darwinian selection. However, it has been remarkably difficult
not only to detect minute selective advantages caused by mo-
lecular changes (1), but also to find genetic systems in which
evolutionary hypotheses can be tested experimentally (2). In the
absence of proper experimental systems, molecular adaptation
in higher eukaryotes has been inferred mostly by using statistical
methods (for examples, see refs. 3–5). For several cases, how-
ever, ancestral molecules have been engineered, allowing studies
of functional changes in the past (6). These analyses demonstrate
that functional changes actually occurred, but they do not
necessarily mean that the new characters were adaptive (7). To
complicate the matter further, evolutionary changes are not
always unidirectional and ancestral phenotypes may reappear
during evolution (8, 9). One effective way of exploring the
mechanisms of molecular adaptation is to engineer ancestral
molecules at various stages of evolution and to recapitulate the
changes in their phenotypes through time. To date, the molec-
ular analyses of the origin and evolution of color vision produced
arguably ‘‘the deepest body of knowledge linking differences in
specific genes to differences in ecology and to the evolution of
species’’ (10). The study of dim-light vision provides another
opportunity to explore the adaptation of vertebrates to different
environments.

Results
Rhodopsins. Dim-light vision in vertebrates is mediated by rho-
dopsins, which consist of a transmembrane protein, opsin, and a
chromophore, 11-cis-retinal (11). By interacting with different
opsins, the identical chromophores in different rhodopsins de-

tect various wavelengths of light (reviewed in ref. 12). To explore
the molecular basis of the spectral tuning in rhodopsins, in vitro
assay-based mutagenesis experiments are necessary, in which the
wavelengths of maximal absorption (�maxs) can be measured in
the dark (dark spectra) and/or by subtracting a spectrum mea-
sured after photobleaching from a spectrum evaluated before
light exposure (difference spectra) (for example, see ref. 13). So
far, the �maxs of contemporary rhodopsins measured by using the
in vitro assay vary between 482 and 505 nm (refs. 12 and 14 and
references therein). By using another method, microspectropho-
tometry (MSP), the rhodopsin of a deep-sea fish, shining loose-
jaw (Aristostomias scintillans), has also been reported to have a
�max of 526 nm (15).

To examine whether these �maxs represent the actual variation
of the �maxs of rhodopsins in nature, we isolated the rhodopsins
of migratory fish [Japanese eel (Anguilla japonica) and its close
relative Japanese conger (Conger myriaster)], deep-sea fish [Pa-
cific blackdragon (Idiacanthus antrostomus), Northern lampfish
(Stenobrachius leucopsarus), shining loosejaw (Aristostomias
scintillans), scabbardfish (Lepidopus fitchi), and Pacific viperfish
(Chauliodus macouni)], and freshwater bluefin killifish (Lucania
goodei), which live in diverse light environments (www.fishbase.
org) [see supporting information (SI) Methods and Fig. S1]. The
eel has two paralogous rhodopsins (EEL-A and -B), as do conger
(CONGER-A and -B) and scabbardfish (SCABBARD-A and
-B), whereas the others use one type of rhodopsins (BLACK-
DRAGON, LAMPFISH, LOOSEJAW, VIPERFISH and BFN
KILLIFISH) (16) (see also SI Result 1).

In the in vitro assay, the �maxs of the dark spectra are more
reliable than those of difference spectra (SI Result 2) and,
therefore, unless otherwise specified, the �maxs refer to the
former values throughout the paper. The �maxs were determined
for EEL-A (500 nm), EEL-B (479 nm), CONGER-A (486 nm),
CONGER-B (485 nm), SCABBARD-A (507 nm), SCAB-
BARD-B (481 nm), and BFN KILLIFISH (504 nm) (SI Result
2). The �maxs for LAMPFISH and VIPERFISH could not be
evaluated, but those of their difference spectra were 492 and 489
nm, respectively. Neither dark spectra nor difference spectra
were obtained for BLACKDRAGON and LOOSEJAW. How-
ever, a mutant pigment, which is modeled after LOOSEJAW,
has a difference spectrum �max of 526 nm (see Molecular Basis
of Spectral Tuning). Hence, the range of �480 –525 nm seems to
represent the �maxs of rhodopsins in vertebrates reasonably well.

Author contributions: S.Y. designed research; S.Y., T.T., H.Z., and L.B. performed research;
S.Y. analyzed data; and S.Y. wrote the paper.

The authors declare no confict of interest.

This article is a PNAS Direct Submission.

Data deposition: The sequences reported in this paper have been deposited in the GenBank
database (accession nos. EU407248 –EU407253).

†To whom correspondence should be addressed at: Department of Biology, Rollins Re-
search Center, Emory University, 1510 Clifton Road, Atlanta, GA 30322. E-mail:
syokoya@emory.edu.

This article contains supporting information online at www.pnas.org/cgi/content/full/
0802426105/DCSupplemental.

© 2008 by The National Academy of Sciences of the USA

13480 –13485 � PNAS � Septmber 9, 2008 � vol. 105 � no. 36 www.pnas.org�cgi�doi�10.1073�pnas.0802426105

D
o
w

n
lo

a
d
e
d
b

y
g
u
e
st

o
n
D

e
ce

m
b
e
r

2
9
,
2
0
2
1

The Ecology of Dim-Light and Deep-Sea Vision. One of the critical
times for the survival of animals in shallow water and on land is
at twilight when the most abundant light falls between 400 and
500 nm (17). Many fish, amphibians, birds, and mammals that
live in these environments use rhodopsins with �maxs of �500 nm
(12). In contrast, in deep water, the distribution of light is much
narrower at �480 nm (18). Mature conger, mature eel, thorny-
head, and coelacanth all live at the depths of 200 –1,800 m
(www.fishbase.org). Our data and that of others (19, 20) show
that these fishes achieve their dim-light vision by using rho-
dopsins with �maxs of �480 nm. Because of their �maxs and
specific light environments, the two groups of rhodopsins may be
classified simply as ‘‘surface’’ and ‘‘deep-sea,’’ respectively.
Despite being active at much deeper depths of 3,000 – 4,000 m
(www.fishbase.org), however, Northern lampfish and Pacific
viperfish use rhodopsins with �maxs of �490 nm. The higher
�maxs can be explained by their upward migration at night, and
LAMPFISH and VIPERFISH can be regarded as ‘‘intermedi-
ate’’ rhodopsins. Then, through the use of far-red (�700 nm)
bioluminescence to create an artificial light environment, shining
loosejaw achieves dim-light vision with the ‘‘red-shifted’’ rho-
dopsins (15).

Japanese eel spawns in the deep sea, the young adults migrate
into freshwater, and the mature fish return to the deep sea for
reproduction. Similarly, Japanese conger spawns in the deep sea,
their larvae hatch near the coast, but they live only in the sea
(www.fishbase.org). For their dim-light vision, young and adult
eels use EEL-A and EEL-B, respectively, whereas congers use
only CONGER-B; CONGER-A is expressed in the pineal
complex (16). The �max of EEL-A (500 nm) ref lects the shallow
freshwater environment, whereas those of EEL-B and CON-
GER-B (480 – 485 nm) match with their light environments in the
deep-sea. The �max of CONGER-A (486 nm) is similar to those
of 470 – 482 nm in the pineal gland-specific pigments of American
chameleon, pigeon, and chicken (reviewed in ref. 12). Hence,
EEL-A is a surface rhodopsin and EEL-B and CONGER-B are
deep-sea rhodopsins. CONGER-A does not ref lect the deep-sea
environment directly; however, because of its �max, CONGER-A
may also be classified as a deep-sea rhodopsin.

Based on considerations of ecology, life history, and �maxs of
rhodopsins, dim-light vision can be classified into biologically
meaningful deep-sea, intermediate, surface, and red-shifted
vision (see SI Result 3 for detailed discussion of the classifications
of other rhodopsins). The corresponding rhodopsins have �maxs
of 479 – 486, 491– 496, 500 –507, and 526 nm, respectively, estab-
lishing the units of possible selection. Consequently, it is possible
that selective force may be able to differentiate even 4 –5 nm of
�max differences of rhodopsins.

Ancestral Rhodopsins. Based on the composite phylogenetic tree
(Fig. 1) (see also SI Result 4 and Fig. S2) of the 11 newly
characterized rhodopsins and 27 others from a wide range of
vertebrate species, we inferred the amino acid sequences of
ancestral rhodopsins. Because most of their rhodopsin genes
have been sequenced partially (Fig. S3), squirrelfish (21) and
cichlid (14) rhodopsins were first excluded from this inference.
The amino acids inferred by using the Jones, Taylor, and
Thornton and Dayhoff models of amino acid replacements of the
PAML program (3) are highly reliable (SI Results 5 and 6, and
Tables S1 and S2).

By introducing a total of 137 amino acid changes into various
rhodopsins (Table S3), we then engineered pigments at nodes
a– k (pigments a– k). The in vitro assays show that pigments a– d
and f– h have �maxs of 501–502 nm, whereas others have �maxs of
496 nm (pigment i) and 482– 486 nm (pigments e, j, and k), which
are also highly reliable (SI Results 5 and 6).

Fig. 1. A composite tree topology of 38 representative rhodopsins in ver-
tebrates. Numbers in ovals are �maxs evaluated from MSP (*), dark spectra, and
difference spectra (†). The numbers in white, blue, black, and red ovals
indicate surface, intermediate, deep-sea, and red-shifted rhodopsins, respec-
tively, whereas those in rectangles indicate the expected values based on the
mutagenesis results. Because of their incomplete data, the amino acid se-
quences of pigments a– k have been inferred by excluding the squirrelfsh,
bluefn killifsh, and cichlid rhodopsins. S-PUN and S-XAN are classifed as
intermediate rhodopsins because of their expected �maxs, but currently avail-
able data are ambiguous. Red- and blue-colored amino acid replacements
indicate the color of the shifts in the �max. The �max of avian ancestral pigment
shows that of the ancestral Archosaur rhodopsin (39). ND, the �max could not
be determined.

Molecular Basis of Spectral Tuning. Because of the interactions
between the 11-cis-retinal and various amino acids, the �max
shifts caused by mutations in the opposite directions are often
nonsymmetrical (22–24). Hence, to understand the evolutionary
mechanism that has generated the various �maxs of rhodopsins in
nature, we must analyze ‘‘forward’’ amino acid replacements that
actually took place in specific lineages. The reconstruction of
multiple ancestral pigments opens an unprecedented opportu-
nity to study the effects of such forward amino acid replacements
on the �max shift in different lineages.

At present, certain amino acid changes at a total of 26 sites are
known to have generated various �maxs of rhodopsins and other
paralogous visual pigments in vertebrates (25). The amino acid
sequences of the 38 representative rhodopsins differ at 11 of the
26 sites (Fig. 2, second column). Among these, the specific amino
acid replacements at 46, 49, 52, 93, 97, 116, and 164 are unlikely
to have been involved in the spectral tuning (SI Result 7). These
and other amino acid site numbers in this paper are standardized

EV
O

LU
TI

O
N

D
o
w

n
lo

a
d
e
d
b

y
g
u
e
st

o
n
D

e
ce

m
b
e
r

2
9
,
2
0
2
1

Yokoyama et al. PNAS � Septmber 9, 2008 � vol. 105 � no. 36 � 13481

D
o
w

n
lo

a
d
e
d
b

y
g
u
e
st

o
n
D

e
ce

m
b
e
r

2
9
,
2
0
2
1

11122 1111223
44589912669 90899581
69233762412 62345397

pigment a FLFDTTFEAFA YYMLKMTM
pigment b ……….. ……..
pigment c .I……… ……..
pigment d .I……… …PN…
pigment e .I……..S ……..
pigment f ……….. …RA…
pigment g ……….. …RA…
pigment h ……….. …RA…
pigment i ……….. VFLRA…
pigment j …N……S …RA…
pigment k …N.S….S …RA…
pigment m ……….. ……..
CONGER-A .I……..S ..LRA…
EEL-A .I……… …PN…
CONGER-B .I……..S ..L…..
EEL-B T..N……S ……..
CAVEFISH …….I.Y. …RA…
GOLDFISH .I……… …RP…
ZEBRAFISH .I……… …RT…
N-SAM .I…SSM.Y. …RA…
N-ARG .I…SSM.Y. …RA…
S-PUN .I…SSM… …RA…
S-MIC .I…..M… …RA…
S-DIA .I…..M… …RA…
S-XAN .I…..M… …RA…
N-AUR .I…S.M..S …RA.A.
S-SPI .I…S.M.YS …RA.A.
S-TIE .I…S.M.YS …RA.A.
M-VIO .I….SMG.. …RA…
M-BER .I….SMGG. …RA…
BFN KILLIFISH .I……… …RA…
O-NIL ……….. …RA…
X-CAU ……….S ..LRA…
THORNYHEAD L..N……S …RA…
LAMPFISH ……MQ… …RA…
SCABBARD-A .I…….Y. VFLRA…
SCABBARD-B …N……S VFLRA…
VIPERFISH …N.S….S …RA.A.
BLACKDRAGON …N.S….S …RA.A.
LOOSEJAW …N…..YI ..FRAL.I
COELACANTH ….VS.Q..S ……G.
CLAWED FROG L.LNV…… …….L
SALAMANDER L..NVS….. ……..
CHAMELEON L..N……. …PT…
PIGEON M………. ……..
CHICKEN M………. ……..
ZEBRA FINCH M………. ……..
BOVINE LM……… …PH…
DOLPHIN LV.N……S …SR…
ELEPHANT LV.N……. ……..

Fig. 2. Amino acids at the 11 previously known (25) and newly found critical
residues of rhodopsins. The numerical column headings specify the amino acid
positions, and the third column describes the newly discovered critical resi-
dues. Shades indicate amino acid replacements that are unlikely to cause any
�max shifts (SI Result 7). Dots indicate the identity of the amino acids with those
of pigment a. The ancestral amino acids that have a posterior probability of
95% or less are underlined.

by those of the bovine rhodopsin (BOVINE). From the mu-
tagenesis results (SI Result 8 and Table S4), four key observations
of evolutionary significance emerge.

First, the �maxs of most contemporary rhodopsins can be
explained largely by a total of 15 amino acid replacements at 12
sites. Namely, significant �max shifts have been caused by 4 of the
11 currently known sites (D83N, E122M, E122Q, F261Y, A292S,
and S292I) as well as newly discovered sites Y96V, Y102F,
E122I, M183F, P194R, N195A, M253L, T289G, and M317I (Fig.
2, third column). Therefore, the functional differentiation of
vertebrate rhodopsin was caused mostly by only �3% of 354
amino acid sites.

Second, 4 of the 15 critical amino acid replacements occurred
multiple times during rhodopsin evolution: D83N (seven times),
A292S (nine times), F261Y (five times), E122Q (two times), and
D83N/A292S (five times) (Fig. 1). Such extensive parallel
changes strongly implicate the importance of these and other
amino acid replacements at the 12 sites in the functional
adaptation of vertebrate dim-light vision.

Third, we uncovered new types of amino acid interactions.
A292S usually decreases the �max of rhodopsin by �10 nm (12)
(see also pigments c and g in Table S4). Much to our surprise,
when A292S was introduced into pigment d (Fig. 1), it did not
decrease the �max at all. However, when the reverse mutation,
S292A, was introduced into CONGER-A, a descendant of
pigment d, the mutant rhodopsin increased the �max by 12 nm,
explaining the �max of pigment d reasonably well. From the latter
analysis alone, we may erroneously conclude that the �max of
CONGER-A was achieved by A292S; instead, it was achieved
purely by the interaction of three amino acid replacements
(P194R, N195A, and A292S) (SI Result 8). Moreover, F261Y in
pigment b increases the �max by 10 nm, and CAVEFISH, a
descendant of pigment b, should have a �max of �510 nm, but the
observed value is 504 nm (12), where the effect of F261Y was
countered by E122I (for more details, see SI Result 8).

To study the molecular basis of spectral tuning, quantum
chemists analyze the interactions between the 11-cis-retinal and
amino acids that are located in the retinal binding pocket, within
4.5 Å of the 11-cis-retinal (26, 27). The residue 292 is �4.5 Å
away from the 11-cis-retinal and is very close to the functionally
critical hydrogen bonded network (28). However, the tertiary
structure of the bovine rhodopsin (28) shows that the residues
194 and 195 are �20 Å away from residue 292 and are not even
in the transmembrane segments. This magnitude of �max shift
and the distance of interacting amino acids are totally unex-
pected.

The fourth significant observation is that the �max shifts of
rhodopsins can be cyclic during vertebrate evolution. In partic-
ular, F261Y reversed the direction of �max shift four times (Fig.
1). If F261Y preceded E122I in the CAVEFISH lineage, then the
functional reversions occurred five times.

Paleontology, Ecology, and Habitats. The ancestors of bony fish most
likely used rhodopsins with �maxs of �500 nm (Fig. 1). What types
of light environment did these ancestors have? The origin of many
early vertebrate ancestors is controversial, but that of bony fish
ancestors is clear (29). The fossil records from late Cambrian and
early Ordovician, �500 Mya, show that the ancestors of bony fish
lived in shallow, near-shore marine environments (30 –32). There-
fore, pigment a must have functioned as a surface rhodopsin and its
�max would be consistent with that role. Interpolating from the
ancestral and contemporary rhodopsins, it is most likely that
pigments b– d and f– h (�max � 501–502 nm) were also surface
rhodopsins, pigment i (496 nm) was an intermediate rhodopsin, and
pigments e, j, and k (480 – 485 nm) were deep-sea rhodopsins (Fig.
1). From their predicted �maxs, it is also likely that pigments q, r, s,
and v were intermediate rhodospins and pigment u was a deep-sea
rhodopsin (Fig. 1).

Based on the four types of dim-light vision, vertebrates show
six different evolutionary paths (Fig. 1). First, surface vision has
been maintained in a wide range of species, from eels to
mammals. Second, the transition of surface 3 intermediate
vision also occurred in a wide range of species. Third, many
deep-sea fish have achieved the directed transitions of surface 3
intermediate 3 deep-sea vision. The three additional changes
are surface 3 intermediate 3 surface vision (some squirrelfish),
surface 3 intermediate 3 deep-sea 3 intermediate vision
(some squirrelfish and Pacific viperfish), and surface 3 inter-
mediate 3 deep-sea 3 red-shifted vision (shining loosejaw), all
showing that the evolution of dim-light vision is reversible.

Molecular Evolution. In vertebrate rhodopsins, several amino acid
replacements occurred multiple times and, furthermore, the
biologically significant �max shifts occurred on at least 18 sepa-
rate occasions (Fig. 1). These observations strongly suggest that
the 15 amino acid changes have undergone positive selection. To
search for positively selected amino acid sites, we applied the

13482 � www.pnas.org�cgi�doi�10.1073�pnas.0802426105 Yokoyama et al.

Table 1. Results from the NEB and BEB analyses

Rhodopsins Model* NEB BEB

Squirrelfsh M2a 50 162 213 214 50 162 213 214
M8 37 50 162 213 214 37 50 112 162

213 214 217
Squirrelfsh and other fsh M2a 162 212 162, 212

M8 162, 212 162, 212
Coelacanth and tetrapods M2a None None

M8 None None
All M2a None None

M8 None None

Sites with P � 0.01 levels are in bold.
*The null models and other parameters are given in Table S5.

naive-empirical-Bayes (NEB) and Bayes-empirical-Bayes (BEB)
approaches of maximum-likelihood-based Bayesian method (3,
33) and parsimony method (4) to four sets of rhodopsin genes:
(i) 11 squirrelfish genes; (ii) the squirrelfish rhodopsins and the
bluefin killifish, cichlid, and deep-sea fish genes, excluding the
coelacanth gene; (iii) the coelacanth and 9 tetrapod genes; and
(iv) all 38 genes (Fig. 1).

Using the parsimony method, we could not find any positively
selected amino acid sites. Using the Bayesian methods, however,
a total of eight positively selected sites (positions 37, 50, 112, 162,
212, 213, 214, and 217) are predicted (Table 1). The Bayesian
results reveal two characteristics. First, the positively selected
amino acid sites are predicted only for relatively closely related
genes, involving squirrelfish genes, but they disappear as more
distantly related genes are considered together. It is also sur-
prising that none of these predicted sites coincide with those
detected by mutagenesis experiments. Second, different amino
acids at these predicted sites do not seem to cause any �max-shifts
(SI Result 9, Tables S5–S8, and Fig. S4).

Considering 17 closely related cichlid rhodopsin genes, 26
positively selected sites have also been predicted (34). Again,
none of the different amino acids at these sites seem to cause any
�max shifts and, furthermore, when we add the other 23 genes
(the eel, conger, cavefish, goldfish, zebrafsh, deep-sea fish, and
tetrapod rhodopsin genes in Fig. 1) in the Bayesian analyses, all
positively selected sites disappear (SI Result 9)! Why can posi-

rhodopsins. According to their �maxs and light environments,
rhodopsins are classified into four groups: deep-sea (�480 – 485
nm), intermediate (�490 – 495 nm), surface (�500 –507 nm), and
red-shifted (�525 nm) rhodopsins. Our mutagenesis results
establish five fundamental features of molecular evolution that
cannot be learned from the standard statistical analyses of
protein sequence data.

First, mutagenesis experiments can offer critical and decisive
tests of whether or not candidate amino acid changes actually
cause any functional changes. Second, the same amino acid
replacements do not always produce the same functional change
but can be affected by the amino acid composition of the
molecule. Therefore, the likelihood of parallel amino acid
replacements, which may or may not result in any functional
change, can overestimate the actual probability of functional
adaptations (SI Result 9).

Third, similar functional changes can be achieved by differ-
ent amino acid replacements; for example, D83N/A 292S,
P194R/N195A/A 292S, and E122Q decrease the �max by 14 –20
nm (Table S4). Thus, by simply looking for parallel replace-
ments of specific amino acids, one can fail to discover other
amino acid replacements that generate the same functional
change, thereby underestimating the chance of finding func-
tional adaptations.

EV
O

LU
TI

O
N

tively selected sites be inferred more often in closely related
genes than in distantly related genes? When nucleotide changes
occur at random, the proportions of nonsynonymous and syn-
onymous mutations are roughly 70% and 30%, respectively.
Hence, under neutral evolution, or even under some purifying
selection, closely related molecules can initially accumulate
more nonsynonymous changes than synonymous changes. How-
ever, as the evolutionary time increases, synonymous mutations
will accumulate more often than nonsynonymous mutations
(35). The differential rates of synonymous and nonsynonymous
nucleotide substitutions during evolution may explain the pre-
diction of false-positives among the relatively closely related
rhodopsins. Or, such inferences may also be affected by the
statistical procedures (36).

As we saw earlier, D83N, Y96V, Y102F, E122I, E122M,
E122Q, P194R, N195A, and A292S decreased the �max, whereas
M183F, M253L, F261Y, T289G, S292I, and M317I increased it.
These changes are located within or near the transmembrane
segments, but most of the amino acid replacements at the

COOH

NH2

D
o
w

n
lo

a
d
e
d
b

y
g
u
e
st

o
n
D

e
ce

m
b
e
r

2
9
,
2
0
2
1

remaining 191 neutral sites are scattered all over the rhodopsin
Fig. 3. Secondary structure of BOVINE (26) with a total of 203 naturallymolecule (Fig. 3).
occurring amino acid replacements in the 38 vertebrate rhodopsins, where
seven transmembrane helices are indicat

Levels of naturally occurring DNA
polymorphism correlate with
recombination rates in
D. melanogaster
David J. Begun & Charles F. Aquadro

Secti on of Genetics and Development, Biotechnol ogy Buildi ng,
Cornell University, Ithaca, New York 14853 – 270 3 . USA

Two genomic regions with unusally low recombination rates in
Drosophila melanogaster have normal levels of divergence but
greatly reduced nucleotide diversity’· 2, apparently resulting from
the fixation of advantageous mutations and the associated hitch­
hiking effect 3.4. Here we show that for 20 gene regions from across
the genome, the amount of nucleotide diversity in natural popula­
tions of D . melanogaster is positively correlated with the regional
rate of recombination. This cannot be explained by va riation in
mutation rates and/or functional constraint, because we observe
no correlation between recombination rates and DNA sequence
divergence between D. melanogast er and its sibling species, D.
simulans. We suggest that the correlation may result from genetic
hitch-hiking associated with the fixation of advantageous mutants .
Hitch-hiking thus seems to occur over a large fraction of the
Drosophila genome and may constitute a major constraint on levels
of genetic variation in nature.

T a ble I summ a riz es le v el s o f DNA v a ri ati o n a nd inte rs pe cifi c
di ve rg enc e ( b et wee n D. m ela nogaste r a nd D. simulan s) wh ere
av ail a bl e. T he se estimate s of D NA pol y morphi sm are d eri ve d
fr o m re stri ction site survey s w ith o ne e x cepti o n ( cu bitu s int errup­
rus) a nd th erefore a r e estimate s o f average l eve ls of va ri ati o n
ov er 13 to 65 kil o ba ses ( kb ) from eac h ge n e re g ion. T o explor e
the rel ati o n ship b et we en levels o f D N A sequen ce v a ri a tion a nd
r eco mbin ati o n ra tes we co mp are d estim at es o f nu cle otid e di ve r ­

5 sit y ( w ) and th e coe ffi c i e nt of ex ch an ge\ a mea sur e of r ec ombi ­
n ation ra t e p er p h ys i ca l d ista n ce .

E
0.014

.?:­
·;;;
ID 0 .010
>
i5
(I)
1::> 0.006
-~
(I)

u
::,
z 0.002

•

O l!H~~~~~~~~-~~~-~~~~~~~~–

0 0.02 o.o• 0.06 o.oa 0 .1
Coefficient of exchange

FIG. 1 Scatte rplot of nucleoti de divers ity ( 7T) versus coef ficient of exchange
in D. melanogaster. Autoso mal and X-linked genes are represented by
hatched and closed circles. respect ively. To make autosoma l and X-linked
genes direc t ly compa rable we made the simplifying assu mpti on of equal
numbers of males and fema les. Then, under neutrality , nucleot ide heterozy­
gos it y is an estimate of 3N µ. fo r X-linked genes and 4Nµ. fo r autosomal
genes (N is the e ff ective populat ion size and µ. is the neut ral mut ation
rate) . Therefo re, befo re doing regress ion or co rrelat ion analyses. we mult i­
plied est imates of 7T from X-linked gene regions by fou r-th irds . The recombi­
nat ion rates est imated by the coefficie nt of exchange are for fema les. An
X-linked gene reg ion spends two -th irds of t he t ime in fema les (where it can
recombine) and only one-thi rd of the time in males (where it canno t recom­
bine ). wherea s an autosome spends half it s t ime in fema les and half in
males. Therefore , we multi plied the coeff icient of exchange fo r autos omal
genes and X-linked regio ns by one-half and two- t hirds. respective ly.
Regressio n line is indicated by a solid line.

NATURE · VOL 356 · 9 APRIL 1992

LETTERS TO NATURE

TABLE 1 Coefficient s of exchange and nucleotide heterozygosities in D. me/anogaster
and divergence with D. simulans

Coefficient
Gene region of exchange ” Divergence Reference

Chromosome I (X)
yello w.achaete (y, ac) 0.0045 0.001 0.054 1
phosphogluconate

dehydrogenase gene
(Pgd) 0.0154 0.003 0.029 1

zest e- tko (z, tko) 0.0222 0.004 12
period (per) 0.0520 0.001 0.050 1
white (w) 0.1400 0.009 13
Notch (N) 0.1212 0.005 14
vermilion (v) 0.0590 0.001 0.047 (D.J.B. and C.F.A.,

unpublished
results)

forked ( f) 0.0455 0.002 15
glucose-6.phosphate

dehydrogenase gene
(Z w) 0.0485 0.001 16

suppre ssor of fork ed
(sul f)) 000 50 0.000 15

Chromosome II
sn-glycerol 3-phosphate

dehydogenase gene
(Gpdh) 0.0800 0.008 17

alcohol dehydrogenase
gene (Adh) 0.0647 0.006 0.045 18,19

DOPA decarboxylase IC.FA et al ..
gene (Ode) 0.0184 0.005 unpublished

results)
amylase gene (Amy) 0.0435 0.008 20
Punch (Pu) 0.0718 0.004 17

Chromosome Ill
esterase-6 gene (Est-6) 0.0604 0.005 21
metallothionein-A gene

(MCnA) 0.0083 0.001 0.072 22
heat.shock protein. 70A

gene (Hsp70A) 0.0069 0.002 0.023 23, 24
rosy (ry) 0.0471 0.003 0.050 25

Chromosome IV
cubitus in terruptus

Dominant (ci°) o• 0.000 0.050 2

Nucleotide diversity (,r) is the average pairwise difference for all pairs of sequences
drawn at random from a population, and can be thought of as heterozygosity per
nucleotide5 . The coefficient of exchange for a gene region was calculated by selecting
two genetically defined loci26 that f lank the region of interest and dividing the distance
in map units between the flanking loci by the number of polytene bands between the
loci6 . The number of polytene bands between loci was determined from Bridge’s maps27_
An important assumpt ion underlying the use of thi s me t ric as an index of recombination
rate is that over large stretches of the genome (for example. 20 to 40 polytene bands),
the average amount of DNA per polytene band is roughly similar between regions.
Available data suggest that this is a reasonable assumption for at least much of the
Drosophila genome 6 28 2 ” .

• The recombination rate on the fourth chromosome is effectively zero30

F i gu re I is a sca tt erpl o t o f w ve r su s th e coeffici ent of e xc h a n ge
fo r 20 ge ne s in D. m ela nogaste r. It i s a p pare n t t h at l e v els o f
nucl eo tid e d ive r si ty i n crease as r ates of r eco mbin a t io n in c r ease.
Varia ti o n in reco mb ina ti on r at es ex p la i ns a l arge fr act ion o f th e
va ri at i on i n nucl eo tid e di v er sity a n d t he null h y p o th es i s th at
th e sl ope i s ze ro i s rej ec t ed with hi g h p ro b a b i l i t y (F , = 16.8,
P = 0.0007). Th e n o n -pa r a m etri c Spea rm a n a nd K e nd all
r egress i o n t es t s are a l so si g nifi ca ntl y diff ere nt from ze ro ( Spea r­
m a n ‘s D = 544, P < 0.0 1; K end a ll ‘s r = 0.437, P < 0 .0 1) . T h e
sa m e co n cl u si o n i s reac h ed w h e n d ata from th e white reg i o n
( w hi ch h as a pa rt i cul a rl y hi g h l eve l o f va ri a t io n in D.
me lanogas ter) is ex clu d ed fr o m the an a l ys i s ( F, = 5.8, P = 0 .03;
Spe ar ma n ‘s D = 544, 0.02 < P < 0.05; K end a ll ‘s r = 0.374, 0.02 <
P < 0.05).

O n e h y p o th es i s to expl a in thi s tr end is tha t gene reg i o n s in
areas of reduced reco mbi na t io n h ave l owe r n eu tra l m ut a tion
rates. Per h aps recom bin a tio n itse l f i s mut age n ic. If th is w er e
t ru e, th en und e r a n e u tra l m o d el th ese gene reg i o n s sho uld a l so
be l ess di verged be t wee n spec i es th an gene reg i ons i n areas of

519
© 1992 Nature Publishing Group

LETTERS TO NATURE

0.08

0.06
QJ
0
C

~ 0.04
Q)
>
o 0.02

e

•

• •

•

0 .1-,-~~~-.–~~ ……. ~—,.–~,-~ ……. ~-.–~–.-

0 0 .010 0.020 0 .030 0.040

Coefficient of exchange

greater recombination rates 7. In Fig. 2 we show a plot of diver­
gence between D. melanogaster and D. simulans versus the
coefficient of exchange, including a ll gene region s for which we
have estimate s of divergence. Clearly, estimate s of divergence
from a larger number of gene regions are desirable. Nevertheless,
the lack of a significant positive regression coefficient (F , =
0.001, P == 0.983) with the available data argues against the
hypothesis that gene regions in areas of low recombination rates
have, on average, lower substitution rates.

Theoretical results show that at the time offixation of a neutral
variant , the amount of linked neutral variation is reduced, and
that the magnitude of the reduction depends on th e recombina­
tion rate 8 • But at a random time (which is an y time a genomic
region in a population is sampled), the average amount of
neutral nucl eotide polymorphism is unaffected by the recombi­
nation rate 8 9• . Therefore, we are unable to arrive at a sa tisfactory
neutral explanation for the patterns seen in Figs 1 and 2.

We propose that the positive correlation between DNA vari­
ation and recombination rate results from th e selective fixation
of advantageous mutants over a significant portion of the
genome. This correlation suggests that levels of neutral va riation
in man y of the gene regions for which variation has been
measured have been reduced by one or more hitch -hiking events.
Provided that a new selectively favoured mutation goes to
fixation before another advantageous mutation arises close to
it, each fixation will be surrounded by a ‘ window ‘ of redu ced
polymorphism, the relative size of which is proportional to the
rate of recombination for that region of the genome 4 • Thus,
where recombin ation rates are very low, ea ch fixation will cause
a wide window of reduced polymorphism, wherea s in reg ion s
of higher recombination, the window will be proportionatel y
smaller . Mo ving along a chromosome tow ards regions of pro ­
gressively lower recombin ation , the windows becom e clos er ,
and ma y begin to overl ap sub stanti ally. Thu s, re gions of low

Received 19 December 1991; accepted 7 Febru.Jry 1992 . 18 . Langley, C. H., Montgomery. E_ & Quat tlebaum, W. F. Proc. natn. Acad. Sci. US.A. 79 , 5631-5635
(1982) .

19 . Aquadro, C. F., Deese, S. F .. Bland, M. M., Langley, C. H. l aurie~Ahlberg, C. C. Genetics U4 , 1 . & Begun, D J. & Aquadr o . C. F. Genetics 129 , 1147 – 1158 (1991 ).
1165-1190 (1986) 2. Berry, A . J .. Ajioka. J. w. & Kr ei tman . M. Genetics 129 . 1111 -11 17 (199 1).

20. Langley. C. H. et al. Genetics 119 , 619 – 629 (1988 ). 3 . Maynar d Smit h. J. & Ha igh, J. Genet . Res . 23, 23 – 35 (197 4).
21. Game. A . Y. & Oakes ho lt. J. G. Genetics 126 , 1021-1031 (1990) . 4 . Kaplan . N. L.. Hudson , R. R. & Langley. C. H. Genetics 123. 887 – 899 (1 989)
22. Lange, B. W .. Langley . C.H . & Step han . W . Genelics 126, 92 1-93 2 (1990) 5. Nei, M . Molecular £vo futionary Genetics (Columbia Univ. Press . 1987)
23 . J. 80. 5350-5354 (1983). 6 . Linds le y , D. L. & Sand ler , L. Phil. Trans. R. Soc B. 277, 295 – 312 (19 77). Leigh Brown. A. Proc. natn. Acad . Sc( US.A.
24 . J. D. 290 , 677 – 682 (198 1). 7. Kimura, M. The Neutral Theory of Molecular Evolution & (Cambridge Univ. Press. 1983 ). Leigh Brown, A. !sh -Horow itz . Natu re
25 . Aq uadro . C. F .. 125, Lado . K. M. & Noon . W A. Genetics 119, 87 5 -888 (19881. 8 . Taji ma, F. Geneti cs 447 -4 54 (1990) .

R. 26 . Ashbumer, M. OrosophHa Genetic Maps (Drosophila Information Service 69. 1991) . 9 . Hudson . R. Theor. Popular . Biol . 23. 183 – 20 1 (1983) .
27. Lindsley. D. L. & Grell, E. H. Gene tic Variations of Drosophila melanogaster (Carnegie Inst itute, 10 Birky. C. W. Jr & Wa lsh . J. B. Proc. natn A cad Sci . US .A 85 , 6414-64 18 (1988)

Washington DC. 196 7) 11. McDonald, J, H. & Kre itman, M. Na ture 351, 65 2- 65 4 (1991) .
28. Sousa. V. Chromosome Maps of Drosophila (CRC, Boca Raton. Florida. 1988) . 12 . Aguade. M ., Miyashita , N. & Langley. C.H. Molec. Biol. E:vol, 6, 123 – 130 (19881 .
29 . J .. & 254 . 221 – 225 (1991) . 13 . Miya sllita . N. & L angl ey. C. H. Genetics 120. Merriam. Ashburner . M.. Hartl , D . L. Ka fatos. F. C. Science 199 – 212 (19881.

w., c. 5. 30 . Hochman. B. in Genetics and Biology of Drosophila Vol. 1b (eds Ashburner, M. & Novitski, E.) 14 . Schaeffer. S. A qua dro , C. F. & Langley, H. Mol ec. Biol. Evol. 30 – 40 (1988) .
903 – 928 (Academ ic. New York, 1976) . 1 5 . Langley, C. H. in Population Biology of Genes and Molecule s (eds Takahata, N. & Crow. J F.)

75 – 91 (Baifuk an. Japan).
16 . Eanes . W. F., Aj ioka, J. W .. Hey, J. & Wes ley, C. Molec. Bio/. £vol. 6, 38 4 – 397 (1989). ACKNOWLEDGEMENTS. We thank M. Nachman f or comments and all members of our laboratory for
17 . Takano . T. S .. Ku sakabe . 5 . & Mukai. T. Genet ics 129, 753 – 761 (1991) . discussion. This work was supported by the NIH and NSF

520

FIG. 2. Scatterp lot of sequence divergence between D. me/anogaster and
D. simulans versus coefficient of exchange in 0. melanogaster. Autosomal
and X-linked genes are represented by hatched and closed circles. respec­
tively . Coefficients of exchange for X-linked and aut osomal regions are
modified as described in Fig. 1 legend. Regression line is indicated by a
solid line.

recombination are ‘hit’ by selective sweep s more often, keeping
polymorphism at a lower average level. Mutations driven to
fixation by meiot ic drive or biased gene conversion would ha ve
similar evo lutionary consequences. Hitch-hiking does not affect
interspecific divergence’° , con sistent with the observed lack of
a corre lation between DNA sequence divergence and recombi­
nation rate .

McDonald and Kreitm a n 11 propo sed that patterns of synony ­
mous and non-synonymous va riation at the Adh locu s in and
between three Drosophila species are incompatible with
neutralit y. They sugg ested that selective fixation of amino -acid
polymorphisms at this locus is the best expl a nation for their
data. Furthermore, they speculate th at selective fixations
occur at a large number of loci. Our results are consistent
with thi s view . However, the number of sele ctively favoured
nucleotides relative to the size of the genome could still be
quite small 4 •

C learly, hitch-hiking and reco mbination rate s do not explain
all of the heterogeneity in levels of var iation across the D.
melanogaster genome. Variation in mutation rate and functional
constraint , as well as several different forms of selection, have
roles in shaping local leve ls of DNA sequence variation . But
the analy sis pre sented here pre se nts the first eviden ce th at hitch ­
hiking driv en b y selective fixation of new mutations may con­
stra in le vels of nucleotide pol ymorphi sm over large portions of
the D. melanoga ster genome. Inferenc e of effective population
size from levels of DNA variation may be co mpromis ed by this
phenomenon.

Much effort is being expended to assemble physical and
genetic maps in sever al species. An unexpected benefit of these
genome mapping projec ts is that it will be possible to examine
wheth e r correla tion s between recombination rates and levels of
DNA variation a re a general phenomenon in natural populations
o f other taxa , including human s. D

NATURE · VOL 356 · 9 APRIL 1992

© 1992 Nature Publishing Group

Research Gate

Proj ec t

Project

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/12069307

Positive and Negative Selection in the DAZ Gene Family

Article in Molecular Biology and Evolution · May 2001

DOI: 10.1093/oxfordjournals.molbev.a003831 · Source: PubMed

CITATIONS READS

74 41

2 authors:

Joseph P Bielawski Ziheng Yang

Dalhousie University University College London

179 PUBLICATIONS 7,211 CITATIONS 310 PUBLICATIONS 70,548 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

An Integrated Microbiome and Genetic Analysis of Pediatric Crohn’s Disease View project

The Origin of Plants: Genomes, Rocks, and Biogeochemical Cycles View project

All content following this page was uploaded by Ziheng Yang on 19 January 2016.

The user has requested enhancement of the downloaded file.

Positive and Negative Selection in the DAZ Gene Family

Joseph P. Bielawski and Ziheng Yang
Department of Biology, Galton Laboratory, University College London, London, England

Because a microdeletion containing the DAZ gene is the most frequently observed deletion in infertile men, the
DAZ gene was considered a strong candidate for the azoospermia factor. A recent evolutionary analysis, however,
suggested that DAZ was free from functional constraints and consequently played little or no role in human sper-
matogenesis. The major evidence for this surprising conclusion is that the nonsynonymous substitution rate is similar
to the synonymous rate and to the rate in introns. In this study, we reexamined the evolution of the DAZ gene
family by using maximum-likelihood methods, which accommodate variable selective pressures among sites or
among branches. The results suggest that DAZ is not free from functional constraints. Most amino acids in DAZ
are under strong selective constraint, while a few sites are under diversifying selection with nonsynonymous/
synonymous rate ratios (dN/dS) well above 1. As a result, the average dN/dS ratio over sites is not a sensible measure
of selective pressure on the protein. Lineage-specifc analysis indicated that human members of this gene family
were evolving by positive Darwinian selection, although the evidence was not strong.

Introduction

Azoospermia is the most common form of infertil-
ity in human males (Shinka and Nakahori 1996). A lo-
cus on the human Y chromosome, the azoospermic fac-
tor (AZF), is believed to contain a gene, or genes, cru-
cial to proper differentiation of male germ cells. The
observation that microdeletions at three different loci of
AZF occur in 5%–15% of infertile men supports this
hypothesis (Ferlin et al. 1999). One of these loci (AZFc)
encodes the Deleted in AZoospermia (DAZ) gene. Be-
cause AZFc is the most frequently observed deletion in
infertile men, it was considered a strong candidate for
the azoospermia factor (Ferlin et al. 1999). Genes from
a number of different pathways, however, are required
for normal spermatogenesis (Elliott and Cooke 1997).

DAZ is located on the Y chromosome, but it is
closely related to the autosomal gene DAZL1. While
DAZL1 is present in all vertebrates, DAZ is found only
in Old World Monkeys. Thus, DAZ [Yq11.23 ] is be-
lieved to have evolved via translocation of DAZL1
[3p24] to the Y chromosome (Saxena et al. 1996; Grom-
oll et al. 1999) some time after the divergence of Old
and New World monkeys; Kumar and Hedges (1998)
dated that divergence to about 40 MYA. After the trans-
location event, DAZ underwent a series of rearrange-
ments and a modifed copy was amplifed, yielding a Y
gene cluster.

DAZ and DAZL1 have a functional role in fertility.
Both DAZ and DAZL1 are expressed exclusively in germ
cells (Cooke et al. 1996; Ruggiu et al. 1997; Gromoll
et al. 1999), and in humans DAZ expression is highest
in spermatogonia (Menke, Mutter, and Page 1997). Ex-
perimental elimination of DAZL1 in mice results in ter-
mination of germ cell development beyond the sper-
matoginial stage (Ruggiu et al. 1997). Moreover, Y-en-
coded human DAZ can compliment the sterile phenotype
of DAZL1 null mice, yielding a partial recovery of sper-

Key words: DAZ, DAZL1, gene family, maximum likelihood,
codon model, positive selection.

Address for correspondence and reprints: J. P. Bielawski, Depart-
ment of Biology, University College London, 4 Stephenson Way, Lon-
don NW1 2HE, United Kingdom. E-mail: j.bielawski@ucl.ac.uk.

Mol. Biol. Evol. 18(4):523–529. 2001
q 2001 by the Society for Molecular Biology and Evolution. ISSN: 0737-4038

matogenesis, which suggests the same or similar target
mRNA for DAZ and DAZL1 during spermatogenesis
(Slee et al. 1999). Although the specifc functions of
DAZ and DAZL1 are unknown, the presence of RNA
recognition motifs suggests that these genes could be
involved in controlling the cell cycle switch from mi-
totic to meiotic cell division (Gromoll et al. 1999); this
cell cycle switch is controlled by RNA-binding proteins
in yeast (Watanabe et al. 1997).

Surprisingly, a recent evolutionary analysis of the
DAZ family (DAZ and DAZL1 genes) indicated a lack
of functional constraints on DAZ. Agulnik et al. (1998)
found a high rate of nonsynonymous substitution, sim-
ilar rates between exons and introns, and similar rates
among the three codon positions. They concluded that
there were no functional constraints on evolution of
DAZ and that patterns of sequence divergence were due
to neutral drift. They hypothesized that Y-linked DAZ
played little role in human spermatogenesis.

The nonsynonymous-to-synonymous rate ratio (dN/
dS) provides a sensitive measure of selective pressure on
the protein. However, when selection pressure varies
among amino acid sites, the average dN/dS ratio might
not be very informative about the evolutionary process-
es affecting the gene. The objective of this study was to
investigate the role of both purifying and positive selec-
tion on the DAZ gene family by using maximum-like-
lihood methods that accommodate differences in selec-
tive pressures among sites (Nielsen and Yang 1998;
Yang et al. 2000). Our fndings indicated that DAZ was
not free of functional constraints and that other expla-
nations for its rapid rate of nonsynonymous evolution
must be considered. There has been considerable debate
as to whether rapid evolution in gene families is caused
by positive Darwinian selection after gene duplication
(Ohta 1993) or by relaxation, but not complete loss, of
functional constraints in redundant genes (Kimura 1983;
Li 1985). In the latter case, a new function might evolve
when formerly neutral substitutions convey a selective
advantage in a novel environment or genetic background
(Dykhuizen and Hartl 1980). We also examined variable
selective pressures among lineages (Yang 1998; Yang
and Nielsen 1998), and our fndings suggest that both

523

D
ow

nloaded from
http://m

be.oxfordjournals.org/ by guest on O
ctober 18, 2015

C (0.001)
(3.47) b

d (1.44)

DAZL1: Mus

a (0.10)

DAZL 1: Macacca

DAZL 1: Human

e (0.35)
r—‘-(o_.3_5_) __ DAZ: Macacca

0.1 L.. __ 9_(_1._1_4) __ DAZ:Human

524 Bielawski and Yang

FIG. 1.—Phylogeny for the DAZ gene family. This topology was
recovered from a maximum-likelihood analysis of nucleotide sequenc-
es and also from least-squares analysis of synonymous divergence.
Branch lengths are proportional to the mean number of nucleotide sub-
stitutions per codon as inferred under model D: free ratios (Yang 1998).
All analyses were conducted using unrooted topologies; this topology
is rooted for convenience. Numbers in parentheses are branch-specifc
v ratios estimated under model D.

models could have played a role in the evolution of the
DAZ gene family.

Materials and Methods
Sequence Data

Two data sets were compiled, refecting a trade-off
between more characters versus more taxa. Data set 1
was composed of 618 bp of DNA sequence (after ex-
clusion of gaps) from fve representatives of the DAZ
gene family (fg. 1). Sequences of DAZL1 were from
Homo sapiens (GenBank accession number U066078),
Macaca mulatta (AF053608), and Mus musculus
(U046694), and sequences of DAZ were from H. sapiens
(NM004081) and Macaca fascicularis (AJ012216). Se-
quences included exons 1–6, A7, C8, and 10. Macaca
fascicularis DAZ (AJ012216) contains an intragenic du-
plication of exons 2–6 and multiple copies of exons 7
and 8. We sampled the 39 copy of exons 2–6, which is
99% similar to the 59 copy. The ‘‘A’’ copy of exon 7
and the ‘‘C’’ copy of exon 8 were sampled because each
predates divergence of the human and Macaca lineages
(Gromoll et al. 1999). Data set 1 was used to investigate
variation in selective pressure among lineages. To study
variation in selective pressure among sites, however,
more sequences were needed. Hence, a second data set
was compiled consisting of 11 members of the DAZ
gene family (fg. 2), but only 291 bp of DNA sequence.
Included in data set 2 were exons 3–5 and portions of
exons 2 and 6. This data set was composed almost en-
tirely of the RNA recognition domain, which spans ex-
ons 2–5. Data set 2 included DAZL1 from Cebua apella
(AF053608), H. sapiens (U066078), M. mulatta
(AF053608), and Papio hamadryas (AF053607); a sin-
gle copy of DAZ from H. sapiens (NM004081), Pan
troglodytes (AF072324), and M. fascicularis
(AJ012216); and two DAZ clones (C1 and C2) from P.
hamadryas (C1: AF07230; C2: AF07321) and M. mu-
latta (C1: AF072322; C2: AF072323). Clones from P.
hamadryas and M. mulatta were divergent copies from
a multicopy DAZ array on the Y chromosome (Agulnick
et al. 1998). Saxena et al. (2000) recently reported that
human DAZ genes occur as a four-copy array in the
AZFc region of the Y chromosome. However, they

found that the four copies differed by only a single,
silent, transition in exon 7A, so only one copy was in-
cluded in our analysis.

Data Analysis

Tree topologies were estimated using maximum
likelihood (ML) under the general time-reversible
(GTR) model with a discrete gamma model (dG) of rate
variation among sites (Yang 1994a, 1994b). Trees also
were estimated by least-squares from synonymous di-
vergences estimated by ML under a codon model of
evolution (Goldman and Yang 1994). The PAUP* com-
puter program (Swofford 2000) was used for conducting
tree searches.

We implemented four nested models of variable se-
lective pressures among branches (Yang 1998; Yang and
Nielsen 1998). Model A was the simplest and assumed
the same v ratio for all branches. Models B and C were
based on the prediction that a gene family evolves under
different selective pressures following gene duplication.
Model B assumed two v ratios: one for the branch pre-
dating the translocation to the Y chromosome (fg. 1;
branch a), and a second for branches postdating the
translocation (branches b–g). Model C assumed three v
ratios: one for branch a, one for DAZL1 branches post-
dating the translocation (branches b–d), and one for all
DAZ branches (branches e–g). Model D (free ratios) as-
sumed an independent v ratio for each branch of a to-
pology and was employed to evaluate the potential for
positive selection in any one branch of the tree.

ML models (Yang and Nielsen 2000; Yang et al.
2000) also permit testing and identifcation of selective
pressures at individual codon sites. We implemented
three such models: M3 (discrete), M7 (beta), and M8
(beta&v). M3 assumed two site classes with the pro-
portions f0 and f1 and ratios v0 and v1 estimated from
the data. M7 assumed that v ratios were distributed
among sites according to a beta distribution. Depending
on parameters p and q, the beta distribution can take a
variety of shapes within the interval (0, 1). M8, an ex-
tension of M7, added an extra class of sites having an
v parameter freely estimated from the data. Positive se-
lection was indicated when an v parameter of M3 or
M8 was .1. The likelihood ratio test was used to com-
pare a one-ratio model (M0) with M3 and to compare
M7 with M8. If there were sites with v . 1, Bayesian
methods were used to calculate the posterior probability
that a site fell into each site class; sites with high prob-
abilities for v . 1 were likely to be under positive Dar-
winian selection (Yang et al. 2000).

All ML analyses of codon models were performed
using the codeml program of the PAML package (Yang
1999). The models employed correction for transition/
transversion rate bias and codon usage bias, features of
DNA sequence evolution that have a signifcant effect
on the estimation of substitution rates (Yang and Nielsen
1998, 2000).

D
ow

nloaded from
http://m

be.oxfordjournals.org/ by guest on O
ctober 18, 2015

A. Maximum likelihood tree from GTR+dr

t

DAZL 1: Cebus ape/la
DAZL 1: Papio hamadryas
DAZL 1: Homo sapiens
DAZL 1: Macacca mulatta

Translocation

0.1

DAZ: Macacca fascicularis
DAZ: Papio hamadryas C2
DAZ: Papio hamadryas C1
DAZ: Macacca mulatta C2
DAZ: Macacca mulatta C1
DAZ: Pan troglodytes

—- DAZ: Homo sapiens

B. Least squares tree from synonymous distances

t

DAZL1: Cebus

DAZL 1: Papio hamadryas
DAZL 1: Homo sapiens
DAZL 1: Macacca mulatta

Translocation

DAZ: Pan troglodytes
DAZ: Homo sapiens
DAZ: Papio hamadryas C1
DAZ: Macacca mulatta C2
DAZ: Macacca mulatta C1

— DAZ: Macacca fascicularis ———– DA Z: Papio hamadryas C2
0.1

Selection in DAZ Gene Family 525

FIG. 2.—Candidate topologies for the DAZ gene family in primates. A, Tree topology recovered from a maximum-likelihood analysis under
the GTR substitution matrix combined with a gamma correction for among-sites rate variation. B, Tree topology recovered from least-squares
analysis of synonymous divergence. Branch lengths are proportional to the mean number of nucleotide substitutions per codon as inferred under
model M8: beta&v (Yang et al. 2000). All analyses were conducted using unrooted topologies; these topologies are rooted for convenience.

D
ow

nloaded from
http://m

be.oxfordjournals.org/ by guest on O
ctober 18, 2015

Results
Variable Selection Pressure Among Lineages

Phylogenetic analyses of data set 1 by different tree
reconstruction methods yielded the same tree topology
(fg. 1), and this topology was used to analyze variable
selective pressures among lineages. Four models (A–D)
were ftted by ML to data set 1 (table 1). The estimate
of v for model A (v 5 0.295) represented an average
over all codon sites and branches. Model A was then
compared with model B, which assumed that selective
constraints changed after translocation of DAZL1 to the
Y chromosome. Twice the difference in their likelihood
scores (2d) was compared with a x2 distribution with
degrees of freedom equal to the difference between
models in number of parameters. This likelihood ratio
test indicated that model B provided a signifcantly bet-
ter ft to these data (2d 5 17.8, df 5 1, P 5 0.000025).

The dN/dS ratio after the translocation, v1 5 0.51, is
signifcantly higher than that prior to the translocation,
v0 5 0.10.

Model B assumed that both autosomal DAZL1 and
the copy that was translocated to the Y chromosome
(DAZ) experienced the same change in selective con-
straints after the translocation event. This simple model
was compared with a more complex model (model C)
in which changes in selective constraints after the trans-
location were allowed to differ between DAZL1 and
DAZ (table 1). Estimates under model C indicate very
different v values for DAZ and DAZL1 after duplication
(table 1). However, the likelihood of model C was not
signifcantly better than that of model B (2d 5 1.12, df
5 1, P 5 0.29).

Because positive selection at any one point in the
phylogeny could have affected our results, we applied

526 Bielawski and Yang

Table 1
Log Likelihood Scores and Parameter Estimates Under Models of Variable Selection Pressures Among Lineages

Model p Parameters for Branches ,

A: One ratio . . . . . . . . . . . . . 1 v0 5 0.295 for all branches 21,442.44
B: Two ratios . . . . . . . . . . . . 2 v0 5 0.102 for branch a 21,433.52

v1 5 0.513 for branches b, c, d, e, f, and g
C: Three ratios . . . . . . . . . . . 3 v0 5 0.103 for branch a 21,432.96

v1 5 0.290 for branches b, c, and d
v2 5 0.574 for branches e, f and g

D: Free ratios . . . . . . . . . . . . 7 v0 5 0.100 for branch a 21,426.40
v1 5 3.474 for branch b
v2 5 0.001 for branch c
v3 5 1.444 for branch d
v4 5 0.350 for branch e
v5 5 0.355 for branch f
v6 5 1.144 for branch g

NOTE.—Analyses were conducted using k as a free parameter and the F61 model of equilibrium codon frequencies. p is the number of branch-specifc v
parameters. v ratios greater than 1 are in bold. Branches are defned in fgure 1.

the free-ratios model (model D) to the same data. The
likelihood score under model D was signifcantly better
than that obtained for model B (2d 5 14.2, df 5 5, P
5 0.014). Branches b, d, and g exhibited v values .1
(table 1). Use of the simpler but less realistic F334
model, which calculates codon frequencies by using
base composition at the three codon positions, produced
similar results. Note that v for branch g was slightly
less than 1 under the F334 model (v6 5 0.99), whereas
it was greater than 1 under the F61 model (v6 5 1.144).

Variable Selection Pressure Among Sites

Phylogenetic analysis of data set 2 under the nu-
cleotide model GTR1dG recovered a topology in which
divergent copies of DAZ from the same species were
not monophyletic, indicating that divergent copies of
DAZ originated in an early amplifcation event and per-
sisted in multiple lineages (fg. 2A). This result is similar
to that obtained in a previous analysis of the DAZ gene
family (Agulnik et al. 1998). We also inferred a tree
topology from synonymous divergences (fg. 2B). This
tree, although different from the tree obtained from the

Table 2

nucleotide analysis, also indicated that some copies of
DAZ originated from an early amplifcation event and
persisted to the present day. Both trees also indicate a
clear bifurcation between all DAZL1 and DAZ sequenc-
es, supporting the hypothesis that a single translocation
event gave rise to the Y-encoded DAZ. To investigate
the impact of tree topology, models of variable v values
among sites were analyzed using both topologies in fg-
ure 2 (table 2). The small size of data set 2 (291 bp; 97
codons) prevented use of the parameter-rich model of
empirical codon frequencies (the F61 model), and the
F334 model was used instead.

The discrete model (M3), which allowed two site
classes with independent v ratios, provided a signifcant
improvement over the one-ratio model (M0) regardless
of the tree topology assumed (table 3). The selective
pressure is not uniform among amino acid sites. Esti-
mates of parameters under M3 suggest that most sites
(95%–97%) are under selective constraint, with v0 5
0.35–0.37, while a few sites (3%–5%) are evolving by
positive selection, with v1 close to 6. Both models, M3
and M8, which allowed for the presence of positively

D
ow

nloaded from
http://m

be.oxfordjournals.org/ by guest on O
ctober 18, 2015

Log Likelihood Scores and Parameter Estimates for Four Models of Variable v’s Among Sites and Two Tree
Topologies

Model Parameter Estimates Positively Selected Sites ,

Tree 1 (fg. 2A)
M0: One ratio . . . . . . . . . . . v 5 0.47 None 2747.19
M3: Discrete . . . . . . . . . . . .

M7: Beta . . . . . . . . . . . . . . .
M8: Beta&v . . . . . . . . . . . .

Tree 2 (fg. 2B)
M0: One ratio . . . . . . . . . . .

v0 5 0.37, f0 5 0.97
v1 5 5.66 (f1 5 0.03)
p 5 0.63, q 5 0.78
p 5 2.2, q 5 3.1, f0 5 0.98
v1 5 12.47 (f1 5 0.02)

v 5 0.47

26, 28, 42

Not allowed
28

None

2742.72

2745.84
2742.70

2758.47
M3: Discrete . . . . . . . . . . . .

M7: Beta . . . . . . . . . . . . . . .
M8: Beta&v . . . . . . . . . . . .

v0 5 0.35, f0 5 0.95
v1 5 5.96 (f1 5 0.05)
p 5 0.30, q 5 0.37
p 5 107, q 5 197, f0 5 0.95
v1 5 5.70 (f1 5 0.05)

26, 28, 42, 91

Not allowed
26, 28, 42, 91

2748.90

2754.99
2748.91

NOTE.— p and q are parameters of the beta distribution. f is the proportion of sites assigned to an individual v category or to a beta distribution with shape
parameters p and q. The proportion f1 (in parentheses) is not a free parameter. Positively selected sites are those with posterior probabilities (P) . 0.50, and those
with P . 0.95 are in bold.

Selection in DAZ Gene Family 527

Table 3
Likelihood Ratio Statistic (2d) for Comparing Models of
Variable v’s Among Sites

M3 vs. M0 M8 vs. M7

Tree 1 (fg. 2A) . . . . . . . . 8.94* 6.82*
Tree 2 (fg 2B) . . . . . . . . 19.14* 12.16*

NOTE.—See table 2 for model parameters.
* Signifcant at the 5% level (x2 5 5.99, df 5 2).5%

selected sites indicated that some variation in selective
pressure was due to positive selection (table 2). Likeli-
hood ratio tests indicated that these models ft the data
better than models in which positively selected sites
were not allowed (table 3). It is also noteworthy that
regardless of model or topology, v values for sites not
subject to positive Darwinian selection were well below
1 (table 2), indicating evolution by purifying selection.

Agulnick et al. (1998) hypothesized that there were
no functional constraints on DAZ sequences. To test this
hypothesis specifcally, we reanalyzed only the DAZ se-
quences of data set 2. The results were consistent with
the previous analysis of data set 2; v values for those
sites not subject to positive Darwinian selection were
well below 1 (e.g., tree 1—discrete model: v0 5 0.43,
f0 5 0.93, v1 5 3.4, f1 5 0.07; beta&v model— p 5
98, q 5 122, f0 5 0.94, v1 5 4.1, f1 5 0.06).

Discussion

Maximum-likelihood analysis of the DAZ gene
family revealed signifcant variation in selective pres-
sures among lineages and among sites. The majority of
sites are clearly subject to purifying selection, with the
nonsynonymous rate being well below the synonymous
rate. A small fraction of sites exhibit nonsynonymous
rates almost six times the synonymous rate, indicating
the action of positive Darwinian selection. Lineage-spe-
cifc analyses indicated that following the translocation
of an autosomal copy of DAZL1 to the Y chromosome,
both loci experienced increased rates of nonsynonymous
substitution. In DAZL1 this was due, at least in part, to
early evolution by positive Darwinian selection. Later,
DAZL1 of M. fascicularis returned to evolution by pu-
rifying selection, whereas DAZL1 of humans continued
to evolve by positive Darwinian selection. Although
there was also an increa

A Single Determinant Dominates the Rate of Yeast Protein Evolution

D. Allan Drummond,* Alpan Raval,�� and Claus O. Wilke§
*Program in Computation and Neural Systems, California Institute of Technology, Pasadena; �Keck Graduate Institute,
Claremont; �School of Mathematical Sciences, Claremont Graduate University; and §Section of Integrative Biology and
Center for Computational Biology and Bioinformatics, University of Texas at Austin

A gene’s rate of sequence evolution is among the most fundamental evolutionary quantities in common use, but what
determines evolutionary rates has remained unclear. Here, we carry out the frst combined analysis of seven predictors
(gene expression level, dispensability, protein abundance, codon adaptation index, gene length, number of protein-protein
interactions, and the gene’s centrality in the interaction network) previously reported to have independent infuences on
protein evolutionary rates. Strikingly, our analysis reveals a single dominant variable linked to the number of translation
events which explains 40-fold more variation in evolutionary rate than any other, suggesting that protein evolutionary rate
has a single major determinant among the seven predictors. The dominant variable explains nearly half the variation in
the rate of synonymous and protein evolution. We show that the two most commonly used methods to disentangle the
determinants of evolutionary rate, partial correlation analysis and ordinary multivariate regression, produce misleading or
spurious results when applied to noisy biological data. We overcome these diffculties by employing principal component
regression, a multivariate regression of evolutionary rate against the principal components of the predictor variables.
Our results support the hypothesis that translational selection governs the rate of synonymous and protein sequence evol-
ution in yeast.

Introduction

A protein’s evolutionary rate, commonly measured by
the number of nonsynonymous substitutions per site in its
encoding gene, is routinely used to characterize functional
importance, detect selection (Nei and Kumar 2000), create
phylogenetic trees (Kurtzman and Robnett 2003), identify
orthologous genes (Wall, Fraser, and Hirsh 2003), and infer
the time of major evolutionary events. However, what
determines a protein’s evolutionary rate has remained the
subject of active speculation and ongoing research (Pál,
Papp, and Hurst 2001; Akashi 2003; Rocha and Danchin
2004).

Recently, studies have found signifcant infuences on
evolutionary rate from many disparate variables: proteins
have been reported to evolve slower if their encoding genes
have a higher expression level in mRNA molecules per cell
(Pál, Papp, and Hurst 2001), if they have a higher codon
adaptation index (CAI) (Rocha and Danchin 2004; Wall
et al. 2005), more protein-protein interactions (higher ‘‘de-
gree’’) (Fraser et al. 2002), shorter length (Marais and Duret
2001), a smaller ftness effect upon gene knockout (higher
‘‘dispensability’’) (Hirsh and Fraser 2001; Yang, Gu, and Li
2003; Zhang and He 2005), or a more central role in inter-
action networks (‘‘betweenness centrality,’’ or simply ‘‘cen-
trality’’) (Hahn and Kern 2005).

Here, we frst demonstrate that the analytical techni-
ques widely used to establish independent roles for many
effects, partial correlation and multiple regression, generate
highly signifcant but entirely spurious effects given noisy
data such as those available for evolutionary analyses.
Then, using a technique which does not suffer from these
problems, we carry out a comprehensive analysis designed
to uncover the major independent correlates of evolutionary
rate in the model eukaryote Saccharomyces cerevisiae. We

Key words: Saccharomyces cerevisiae, evolutionary rate, gene
expression, protein-protein interactions, dispensability, translational
selection.

E-mail: cwilke@mail.utexas.edu.

Mol. Biol. Evol. 23(2):327–337. 2006
doi:10.1093/molbev/msj038
Advance Access publication October 19, 2005

� The Author 2005. Published by Oxford University Press on behalf of
the Society for Molecular Biology and Evolution. All rights reserved.
For permissions, please e-mail: journals.permissions@oxfordjournals.org

determine the number of such correlates, their strength,
and their relationship to the biological variables used in
previous studies. Finally, we ask what these correlates re-
veal about the biological constraints on protein sequence
evolution.

Materials and Methods
Genomic Data

We obtained CAIs and evolutionary rates (nonsynon-
ymous substitutions per site dN, synonymous substitutions
per site dS, adjusted synonymous substitutions dS# [Hirsh,
Fraser, and Wall 2005], and ratios dN/dS and dN/dS#) from
four-way yeast species alignments for 3,036 S. cerevisiae
genes (Wall et al. 2005, supporting information, Table
4). Deletion-strain growth rate data were downloaded from
http://chemogenomics.stanford.edu/supplements/01yfh/fles/
orfgenedata.txt; the average growth rates of the homozy-
gous deletion strains were used as dispensability measure-
ments in our analysis. The fltered yeast interactome data set
(Han et al. 2004) provided interaction network hub types for
199 genes and the number of interactions for 1,379 yeast
genes. The latter data set was used to compute betweenness-
centrality values, which quantify the frequency with which
a network node lies on the shortest path between other
nodes, as described by Hahn and Kern (2005). Genomic
data for Saccharomyces paradoxus and Kluyveromyces
waltii were obtained exactly as described by Drummond
et al. (2005). Genome sequences for Escherichia coli
K12 and Salmonella typhimurium LT2 were obtained from
the Institute for Genomic Research (Peterson et al. 2001),
with orthologs identifed and evolutionary rates computed
exactly as described (Drummond et al. 2005). Gene ex-
pression levels for E. coli measured in mRNAs per cell
in Luria-Bertani (LB) and M9 media were obtained from
Bernstein et al. (2002).

Statistical Analysis

We used R (Ihaka and Gentleman 1996) for statistical
analyses and plotting. The package �pls� was used to perform

D
ow

nloaded from
https://academ

ic.oup.com
/m

be/article/23/2/327/1118974 by guest on 16 D
ecem

ber 2021

328 Drummond et al.

principal component regression. We log transformed all
variables except dispensability. We decided whether or
not to log transform a variable based on whether log trans-
formation led to a higher R2. For those variables that
contained zeros, we added a small constant before the
log transformation, as previously suggested (Wall et al.
2005). This constant was 0.001 for dN, dS#, and dN/dS#
and 10

�7
for betweenness centrality. We scaled the predic-

tor variables to zero mean and unit variance before carrying
out the principal component analysis. In all regression anal-
yses (both against the original predictors and against the
principal components), we determined the statistical signif-
icance levels by starting with the full model and succes-
sively dropping the least signifcant predictor until only
signifcant predictors (P , 0.01) remained.

Results
Correlation and Partial Correlation Analysis

We used the yeast S. cerevisiae to examine the deter-
minants of evolutionary rate because it has been the subject
of many previous analyses (e.g., Pál, Papp, and Hurst 2001;
Fraser 2005) and has an enormous amount of available ge-
nomic, proteomic, and functional data. We frst examined
the raw correlation of six previously assessed biological
variables (expression, CAI, length, dispensability, degree,
and centrality) with protein evolutionary rate, as measured
by the number of nonsynonymous substitutions per site in
the underlying gene. A seventh variable, the number of pro-
tein molecules per cell (‘‘abundance’’), was also considered.
Table 1 shows that all variables except centrality correlated
signifcantly with evolutionary rate, as previously reported.

Expression level strongly correlates with evolutionary
rate, and higher expressed genes have higher CAIs (Akashi
2001), are less dispensable (Gu et al. 2003), more abundant
(Ghaemmaghami et al. 2003), and more likely to be found
in protein-protein interaction experiments (Bloom and
Adami 2003) than lower expressed genes. No inverse rela-
tionships have been posited by which these variables alter
the expression level. Thus, it is imperative to establish
whether these variables play a role independent of expres-
sion level. Following previous analyses (Pál, Papp, and
Hurst 2003; Lemos et al. 2005; Wall et al. 2005), we com-
puted the partial correlation of our seven variables with

Table 1
Partial Correlation Analysis of Seven Putative
Determinants of Evolutionary Rate

Correlation Partial Correlation
Variable X rX,dN rX,dNjgene expression VIF

Gene expression �0.537*** 0 2.72
CAI �0.565*** �0.338*** 2.46
Protein abundance �0.478*** �0.232*** 2.05
Gene length 0.136*** 0.010 1.25
Gene dispensability 0.265*** 0.183*** 1.08
Degree (number of �0.246*** �0.127* 1.70

protein-protein interactions)
Protein centrality �0.098# �0.082 1.64

(frequency on node-node
shortest paths)

#
NOTE.— P , 0.01; * P , 10�3; *** P , 10�9 .

evolutionary rate, controlling for expression level. Table 1
shows that CAI, dispensability, and degree all showed re-
duced but highly signifcant partial correlations, consistent
with previous studies (Hirsh and Fraser 2003; Wall et al.
2005), as did abundance.

Partial Correlations and Noisy Data

What can we conclude from highly signifcant partial
correlations? Yeast expression-level measurements from
multiple groups, even two using the same commercial ol-
igonucleotide array, correlated with coeffcients of only
0.39–0.68 (Coghlan and Wolfe 2000), demonstrating that
expression-level measurements are inaccurate and/or
simply refect the variability of gene expression across
growth conditions and strains. We refer to all such variabil-
ity as noise, regardless of its source. Noisy data are the rule
in genome-wide molecular studies, leading us to explore
what effect noise has on partial correlation analyses. As
a concrete example, CAI is so tightly bound to expression
level that a recent analysis used CAI as its preferred expres-
sion-level measurement (Wall et al. 2005). Might CAI’s
signifcant partial correlation only refect our inability to
control for the true (i.e., evolutionarily relevant) underlying
expression level? More generally, we can ask: what is the
expected partial correlation of two variables, controlling for
a third, when (1) the two variables relate only through de-
pendence on the third ‘‘master’’ variable and (2) all meas-
urements contain noise?

Given these conditions, we derive explicit formulas
for the expected partial correlation, its statistical signif-
cance, and its behavior under various limiting cases in
the Appendix. The expected partial correlation is, in gen-
eral, larger than zero because the full correlation refects
the true underlying master variable’s infuence, while par-
tial correlations can only remove the portion of this infu-
ence that is visible through a noisy measurement (box 1).
We show that, surprisingly, if measurements of an underly-
ing causal variable (e.g., expression level) are noisy, highly
signifcant partial correlations of virtually any strength be-
tween the dependent predictors can be obtained.

As a case in point, dispensability’s role has been vig-
orously debated (Hirsh and Fraser 2003; Pál, Papp, and
Hurst 2003; Wall et al. 2005) with correlation and partial
correlations acting as key analytical tools. Given a model
in which expression level X and noise completely determine
dispensability D and evolutionary rate K (see box 1), what
is the observed partial correlation rDKjX# if we ft variables to
approximately match the observed correlations between X#,
D, and K? As a concrete example, previous reports show
that, using parametric Pearson’s correlations, rX#K ’
�0.6 (Pál, Papp, and Hurst 2001; Wall et al. 2005),
rDK ’ 0.25 (Wall et al. 2005), rDX# ’ 0.2 (Pál, Papp, and
Hurst 2003), and rDKjX# ’ 0.24 (Wall et al. 2005). We can
obtain roughly the reported full correlations and rDKjX# ’

10
�9

0.23 6 0.02, P with 3,000 observations if the true
expression level X is normally distributed with mean 0.5
and standard deviation (SD) 0.25, and the observable
predictors X#, D, and K are equal to X plus zero mean nor-
mally distributed noise with SDs of 0.3, 0.7 and 0.1, respec-
tively. This highly signifcant partial correlation is entirely

D
ow

nloaded from
https://academ

ic.oup.com
/m

be/article/23/2/327/1118974 by guest on 16 D
ecem

ber 2021

A Single Determinant Dominates the Rate of Yeast Protein Evolution 329

Table 2
Results of Principal Component Regression Analysis on Seven Predictors and Five Measures
of Evolutionary Rate for 568 Saccharomyces cerevisiae genes

Principal Components

1 2 3 4 5 6 7 All

Percent variance explained in
dN 42.76*** 0.05 0.50 0.19 0.14 0.47 0.48 44.60***
dS 50.77*** 2.13** 0.88* 0.08 6.55*** 0.37 1.14* 61.92***
dN/dS 24.82*** 0.05 0.67 0.82 0.00 0.00 0.05 26.42***
dS# 6.70*** 0.19 7.31*** 0.26 0.06 0.14 1.25# 15.92***
dN/dS# 42.34*** 0.09 0.13 0.28 0.13 0.40 0.70# 44.07***

Percent contributions
Expression 32.8 1.2 1.5 0.1 1.2 11.2 52.1
CAI 28.3 3.1 8.4 0.9 2.7 17.6 39.0
Abundance 29.2 2.0 1.6 0.3 15.4 51.4 0.1
Length 2.0 1.1 86.4 0.0 2.1 0.3 8.2
Dispensability 1.8 13.0 0.0 84.0 0.0 0.9 0.3
Degree 5.0 36.7 1.9 6.2 38.9 10.9 0.4
Centrality 0.9 42.9 0.3 8.5 39.6 7.7 0.0

# ** ***
NOTE.— P , 0.01; * P , 10�3; P , 10�6; P , 10�9. Bold indicates that the indicated predictor contributes at least 20% to

the indicated component.

spurious: in this model, expression level and random noise
completely determine dispensability. Thus, the observed
statistical relationship between dispensability and evolu-
tionary rate, established by correlation and partial correla-
tion, would arise even if no actual relationship existed
except mutual dependence on noisily measured expres-
sion level.

Multivariate Regression Analysis

Because partial correlation analysis is not applicable to
the problem at hand, what other methods can we use to de-
termine the relative infuence of different predictors on the
rate of evolution? One obvious choice is multivariate re-
gression analysis, a method with the added beneft that
we can look at the infuence of all potential predictor var-
iables at the same time and can eliminate step by step those
predictors that contribute the least to the regression model.
Indeed, several authors have followed this route (Rocha and
Danchin 2004; Agrafoti et al. 2005). Regressing dN simul-
taneously against the seven predictors we consider here, we
fnd that all but centrality make a signifcant contribution to
the regression and that the overall R2 5 0.45.

Unfortunately, ordinary multivariate regression is not
appropriate to analyze the infuence of the various predic-
tors on the evolutionary rate either (box 1). The problem is
that the predictors intercorrelate, while multivariate regres-
sion implicitly assumes that the predictors are statistically
independent. This problem is widely discussed in the sta-
tistical literature, mostly in the context of ‘‘collinear’’ or
‘‘nearly collinear’’ predictors (Gunst and Mason 1977a,
1977b; Mandel 1982; Næs and Martens 1988). The vari-
ance infation factor (VIF) may be used to quantify the de-
gree of predictor collinearity, and table 1 reports VIFs for
our data. These VIFs indicate some collinearity but are not
high enough to raise signifcant concerns. However, for our
toy model (box 1) in which the two predictors refect the
same underlying variable plus noise, the VIFs are only
1.21 in both cases, yet the analysis demonstrates that mul-
tivariate regression and partial correlation break down any-

way. Collinearity and noise work together to undermine
these techniques.

Principal Component Regression Analysis

An alternative approach is to frst identify the indepen-
dent sources of variation in the data, and then determine the
contribution of each biological predictor to each source.
The technique of principal component regression offers
a standard way to carry out such an analysis.

In principal component regression (Mandel 1982),
multiple linear predictors (e.g., expression level and dis-
pensability) are scaled to zero mean and unit variance,
inserted in a matrix, and rotated such that the new coordi-
nate axes point in the directions of greatest predictor var-
iation. The new axes defne variables, called principal
components, which are linear combinations of the original
predictors. Subsequent linear regression of the response
(e.g., dN) on the rotated predictor data yields several pieces
of information per principal component: the proportion of
the response’s variance, R2, explained by the component,
the signifcance of this R2, and the fractional contribution
of each original predictor to the component. Because all
principal components are orthogonal and independent,
the total proportion of response variance explained by
the data is the sum of the component R2 values. Principal
component regression thus circumvents the debilitating
problems of partial correlation and multivariate regression
analyses (box 1) while yielding results which are, in some
ways, easier to interpret.

We carried out principal component regression on the
seven predictors analyzed above. Because the determina-
tion of principal components involves only the predictors
and not the response (i.e., dN or dS), there is only one
set of components and contributions from biological predic-
tors. The regression analysis generates response-specifc
results, in particular, the proportions of variance in dN,
dS, and so on, which each component explains. Table 2
shows numerical data from the analysis of dN and dS using
the seven predictors of expression, CAI, abundance, length,

D
ow

nloaded from
https://academ

ic.oup.com
/m

be/article/23/2/327/1118974 by guest on 16 D
ecem

ber 2021

a * b a * b so * 40 * 40
Expr ssion

30

§ ____
Expr ssion 30 N’ 40 e; N’

-u 30 20
e; 20

(I) -u
C: (I)

·;;; C: 30 10
C. CAI

10 ·;;; CAI § * * * * X 20 C. 0 —-== (I) X
(I) 0 (I) 20 2 3 4 5 u 2 3 4 5 (I)
C: u
ro C:

·;:: ro
ro 10 Abu dance

·;::
> ro 10 >
cf!. cf!. * * * * Ed * 0 — = — — — – 0 ~ — = = E==3

2 3 4 5 6 7 2 3 4 5 6 7

Principal component Principal component

330 Drummond et al.

FIG. 1.—Principal components regression on the rate of protein evo-
lution (dN) in 568 yeast genes reveals a single dominant underlying com-
ponent. (a) Of the seven principal components only one (starred) explained
a statistically signifcant proportion of the variation in dN. This component
explained 43% of the variance, while no other component explained more
than 1%. Expression level, CAI, and protein abundance determined most
of this dominant component (labeled), while the remaining predictors (in
order from top to bottom: length, dispensability, degree, and centrality)
determined ,10% of the component’s R2. See table 2 for numerical data.
(b) A larger data set (1,939 genes) excluding protein-protein interaction
predictors showed the same patterns as in (a).

dispensability, degree, and centrality; fgures 1a and 2a
show these data graphically.

Strikingly, for the rate of protein evolution, dN, one
principal component explained 43% of the variance with
high signifcance, while all other components explained
less than 1% (fg. 1). The single dominant component
was almost entirely (.90%) determined by roughly equal
contributions from three predictors: expression level, abun-
dance, and CAI.

While the causes of dNs variation have remained un-
clear, dS is constrained by translational selection: selection
for preferred codons, which correspond to abundant tRNAs
and are translated faster and more accurately (Akashi 1994,
2001), makes many synonymous changes unfavorable and
thus reduces dS (Hirsh, Fraser, and Wall 2005). Figure 2
shows that the dS results mirror those using dN: the frst
component, which is determined almost entirely by expres-
sion, abundance and CAI, is overwhelmingly dominant
(50.8% of dS variation). A second highly signifcant com-
ponent of modest size (6.6% of dS variation) appears but is
88.4% determined by abundance and CAI. Astonishingly,
the seven biological predictors explain a cumulative 61.9%
of the total variance in dS, with three predictors (expression,
abundance, and CAI) contributing roughly equal amounts
and accounting for 87% of the total variance explained. Be-
cause synonymous sites are thought to be under relatively
weak selection, we would expect random fuctuations
(noise) to contribute a large proportion of variation in
dS, yet our analysis suggests that selective pressures, even
those revealed using noisy data, account for almost two-
thirds of the dS variation among these genes.

The size of the seven-component data set (568 genes)
was severely limited by the requirement for genes having
measures for all seven predictors. In particular, we used
high-quality interactions measurements (Han et al. 2004)
for degree and betweenness centrality; eliminating these

FIG. 2.—Principal components regression on the rate of synonymous
site evolution (dS) in 568 yeast genes reveals a single dominant underlying
component. (a) Seven-predictor variables (see text) yielded seven principal
components, of which six (starred) explained a statistically signifcant pro-
portion of the variation in dS. The dominant component explained 51%
of the variance, while no other component explained more than 7%. See
fgure 1 caption for the breakdown of predictor contributions. (b) A larger
data set (1,939 genes) excluding protein-protein interaction predictors
showed the same patterns as in (a).

measurements, which apparently contribute negligible
amounts to evolutionary rate, more than triples the data
set size to 1,939 genes. We performed the same analysis
on this expanded set and obtained similar results (table 3,
and fgs. 1b and 2b).

To examine the possible effects of assuming a linear
model in the regression, we repeated our analyses using
only data ranks for the predictors and each response.
The results of this nonparametric analysis were virtually un-
changed from the parametric case (data not shown), indi-
cating that little information is contained in the relative
magnitudes of the variables.

It is common practice to interpret dS as the rate of se-
lectively neutral divergence and the ratio dN/dS as the de-
viation of protein evolutionary rate from neutral, putatively
allowing detection of purifying selection or adaptive evo-
lution. We analyzed dN/dS and found trends that were sim-
ilar to those observed in dN and dS alone (tables 2 and 3).
The dominant principal component explained only half the
variation in dN/dS compared to dN or dS, but the reason
seems obvious in light of our results: dN and dS appear
to refect the same underlying selective force, so dividing
one by the other removes much of the shared infuence.
In yeast, as in many other organisms, dS does not refect
neutral divergence but rather divergence constrained by
translational selection for preferred codons, as previous
authors have noted (Hirsh, Fraser, and Wall 2005). These
authors proposed an adjusted measure of dS, denoted dS#,
from which the infuence of codon preference has been
extracted (Hirsh, Fraser, and Wall 2005). We thus analyzed
dS# and dN/dS# (tables 2 and 3), and found that for dS# the
dominance of the frst principal component was obliterated.
While two components (component 1, mostly CAI, expres-
sion and abundance; and component 2, mostly dispensabil-
ity) appeared to make small but possibly meaningful
contributions (R2 . 6%) in the smaller seven-predictor data
set, these contributions were effectively eliminated in the

D
ow

nloaded from
https://academ

ic.oup.com
/m

be/article/23/2/327/1118974 by guest on 16 D
ecem

ber 2021

A Single Determinant Dominates the Rate of Yeast Protein Evolution 331

Table 3
Results of Principal Component Regression Analysis on Five Predictors and
Five Measures of Evolutionary Rate for 1,939 Saccharomyces cerevisiae genes

Principal Components

1 2 3 4 5 All

Percent variance explained in
dN 36.94*** 0.05 0.03 0.22 0.60*** 37.85***
dS 39.33*** 0.73** 0.09 1.93*** 1.92*** 44.01***
dN/dS 22.39*** 0.28 0.21 0.00 0.21 23.10***
dS# 1.26** 2.52*** 2.58*** 0.00 1.54** 7.91***
dN/dS# 37.61*** 0.28 0.00 0.14 1.16** 39.20***

Percent contributions
Expression 33.2 1.7 0.1 24.2 40.8
CAI 31.4 1.0 9.4 9.0 49.2
Abundance 31.3 0.6 0.4 65.8 1.9
Length 2.0 61.0 29.6 0.4 7.0
Dispensability 2.1 35.7 60.5 0.6 1.1

# ** ***
NOTE.— P , 0.01; P , 10�6; P , 10�9. Bold indicates that the indicated predictor contributes at least 20% to the

indicated component.

larger fve-predictor data set (R2 , 3%), even though the
major contributing predictors were still present. This sam-
ple size dependence suggests that the contributions of com-
ponents 1 and 2 are artifacts. Overall, our results are
consistent with the previous claim that dS# has been purged
of the infuence of selection on synonymous sites (Hirsh,
Fraser, and Wall 2005). As additional support, the dN/
dS# regression was nearly instinguishable from that of
dN (tables 2 and 3).

To assess the importance of phylogenetic distance on
our results, we carried out principal component regression
on dN and dS values calculated using two relatives of
S. cerevisiae, S. paradoxus and K. waltii, which diverged
roughly 5 and 100 MYA, respectively (Drummond et al.
2005).

For S. paradoxus, we obtained almost identical results
for dN as for the data of Hirsh, Fraser, and Wall (2005).
However, dS showed a much weaker, though still dominant,
frst component that explained 15% of the dS variance in-
cluding interaction data and 6% without these data, fvefold
more than any other variable. We traced the weaker dS signal
to differences in gene fltering (the smaller data set of Hirsh,
Fraser, and Wall (2005) omits sequences whose gene-level
phylogeny did not match the species-level pattern and se-
quences containing introns and potential frameshifts) and
in codon frequency estimates. Controlling for gene fltering,
the nine-free-parameter codon frequency model used by
Hirsh, Fraser, and Wall (2005) produced a larger signal than
the sixty-free-parameter model used by Drummond et al.
(2005), indicating that analyses of dS may be sensitive to
estimation methodologies (data not shown).

For the distant relative K. waltii, we again obtained
nearly identical results for dN. For the 2,412 genes without
(and 752 genes with) interaction data, one principal com-
ponent determined by CAI, abundance, and expression ex-
plained 41% of the variance in dN, while all other
components explained ,2%. For dS, no dominant compo-
nent emerged, and the best component (mostly expression
and CAI) explained 1.7% of the variance. The lack of any
predictive signal for dS is not surprising because the dS val-
ues relative to K. waltii average more than 14 substitutions

per synonymous site, far beyond the range of reliable esti-
mation. These high dS values may result from a combina-
tion of the large amount of time separating the species,
changes in synonymous pressures, and diffculties in ortho-
log identifcation and alignment. The robust dN results lend
weight to the frst two explanations. We expect that as even
more distant relatives are analyzed, the dN results will be
attenuated by noise, alignment degradation, and phenotypic
changes that must, in some cases, be linked to changes in
relative gene expression levels.

To assess whether the trends we identifed for yeast
extend to other species, we examined evolutionary rates
in 2,605 E. coli genes relative to S. typhimurium. Lacking
global protein abundance, interaction, and dispensability
data for E. coli, we used length, two measures of expression
level, refecting growth in minimal M9 and rich LB media,
and two measures of codon optimization, CAI and the fre-
quency of optimal codons Fop (Ikemura 1985), as predic-
tors. Again, a dominant component emerged which
explained 36% of the dN variance (16-fold more than
any other) and 25% of the dS variance (38-fold more than
any other). Because most of the included predictors are
translation oriented in some way, our results offer no con-
clusion as to the possible infuence of other predictors in
E. coli. However, the remarkable similarity to the yeast re-
sults, including the large portion of variance explained, sug-
gests that similar selective forces have shaped evolutionary
rates in this prokaryotic organism.

Analysis of Binary Variables Using Analysis of
Covariance

In all the above analyses, we found that protein-protein
interactions and gene dispensability showed little or no ap-
parent infuence on the rate of protein evolution (dN) and
synonymous site evolution (dS), contrary to previous
reports (Hirsh and Fraser 2001; Fraser et al. 2002; Fraser
and Hirsh 2004; Wall et al. 2005). Perhaps these measures,
as continuous predictors de

Science
~MAS

D
JSTOR

Four Evolutionary Strata on the Human X Chromosome

Author(s): Bruce T. Lahn and David C. Page

Source: Science , Oct. 29, 1999, New Series, Vol. 286, No. 5441 (Oct. 29, 1999), pp. 964-967

Published by: American Association for the Advancement of Science

Stable URL: https://www.jstor.org/stable/2899501

REFERENCES
Linked references are available on JSTOR for this article:
https://www.jstor.org/stable/2899501?seq=1&cid=pdf-
reference#references_tab_contents
You may need to log in to JSTOR to access the linked references.

JSTOR is a not-for-profit service that helps scholars, researchers, and students discover, use, and build upon a wide
range of content in a trusted digital archive. We use information technology and tools to increase productivity and
facilitate new forms of scholarship. For more information about JSTOR, please contact support@jstor.org.

Your use of the JSTOR archive indicates your acceptance of the Terms & Conditions of Use, available at

Terms and Conditions of Use

American Association for the Advancement of Science is collaborating with JSTOR to digitize,
preserve and extend access to Science

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:10:40 UTC�������������

All use subject to https://about.jstor.org/terms

964

mann, J. R. Ecker, Ce// 72, 427 (1993)] (23). The
largest of 16 independent NPH3 cDNAs was se­
quenced (24) completely (GenBank accession num­
ber AF180390).

11. GenBank searches were accomplished with the
gapped BLAST program [S. F. Altschul et al., Nucleic
Acid Res. 25, 3389 (1997)].

12. The data are available at www.sciencemag.org/
feature/data/ 10423 58.shl

13. Single-letter abbreviations for the amino acid resi­
dues are as follows: A, Ala; C, Cys; D, Asp; E, Glu; F,
Phe; G, Gly; H, His; I, Ile; K, Lys; L, Leu; M, Met; N, Asn;
P, Pro; Q, Gln; R, Arg; S, Ser; T, Thr; V, Val; W, Trp; and
Y, Tyr.

14. T. Patschinsky, T. Hunter, F. S. Esch, J. A. Cooper, B. M.
Sefton, Proc. Natl. Acad. Sci. U.S.A. 79, 973 (1982).

15. The BTB/POZ domain was identified with SMART [J.
Schultz, F. Milpetz, P. Bork, C. P. Ponting, Proc. Natl.
Acad. Sci. U.S.A. 95, 5857 (1998)]. The coiled-coil
structure was identified with COILS [A. Lupas, M. Van
Dyke, J. Stock, Science 252, 1162 (1991)].

16. 0. Albagli, P. Dhordain, C. DeWeindt, G. LeCocq, D.
LePince, Cell Growth Differ. 6, 1193 (1995); L. Ara­
vind and E. V. Koonin,J. Mo/. Biol. 285, 1353 (1999).

17. C. Cohen and D. A. D. Parry, Proteins 7, 1 (1990); A.
Lupas, Trends Biochem. Sci. 21, 375 (1996).

18. Structural analyses were performed with the Protean
program (DNASTAR, Madison, WI).

19. S. Fields, Methods 5, 116 (1993); S. Fields and R.
Sternglanz, Trends Genet. 10, 286 (1994).

REPORTS

20. M. Nagao and K. Tanaka, J. Biol. Chem. 267, 17925
(1992).

21. M. C. Faux and J. D. Scott, Cell 85, 9 (1996); T.
Pawson and J. D. Scott, Science 278, 2075 (1997);
E. A. Elion, Science 281, 1625 (1998).

22. S. D. Choi, R. Creelman, J. Mullet, R. A. Wing, Weeds
World 2, 17 (1995), http://genome-www.stanford.
edu/ Arabidopsis/ww/home.html.

23. J. Sambrook, E. F. Fritsch, T. Maniatis, Molecular Clon­
ing: A laboratory Manual (Cold Spring Harbor Labo­
ratory Press, Plainview, NY, 1989).

24. Sequencing templates were prepared by polymerase
chain reaction and sequenced with an ABl377 auto­
mated sequencer (Perkin-Elmer, Norwalk, CT).

25. Phototropism and hypocotyl growth was assayed as
described previously [E. L. Stowe-Evans, R. M. Harper,
A. V. Motchoulski, E. Liscum, Plant Physio/. 118, 1265
(1998)].

26. E. Liscum and R. P. Hangarter, Plant Cell 3, 685
(1991).

27. J. W. Reed, P. Nagpal, D. S. Poole, M. Furuya, J. Chory,
Plant Cell 5, 147 (1993).

28. C. Bell and J. R. Ecker, Genomics 19, 137 (1994).
29. H.-G. Nam et at., Plant Cell 1, 699 (1989).
30. Information about markers AM40 and AM80 is available

at http:/ /www.biosci.missouri.edu/liscum/newmarkers.
html

31. Y. Nakamura et al., DNA Res. 4, 401 (1997); http://
www.kazusa.or.jp/arabi/chr5/map/24-26Mb.html.

32. Soluble and total microsomal membrane fractions

Four Evolutionary Strata on the
Human X Chromosome

Bruce T. Lahn* and David C. Paget

Human sex chromosomes evolved from autosomes. Nineteen ancestral auto­
somal genes persist as differentiated homologs on the X and Y chromosomes.
The ages of individual X-Y gene pairs (measured by nucleotide divergence) and
the locations of their X members on the X chromosome were found to be highly
correlated. Age decreased in stepwise fashion from the distal long arm to the
distal short arm in at least four “evolutionary strata.” Human sex chromosome
evolution was probably punctuated by at least four events, each suppressing
X-Y recqmbination in one stratum, without disturbing gene order on the X
chromosome. The first event, which marked the beginnings of X-Y differenti­
ation, occurred about 240 to 320 million years ago, shortly after divergence of
the mammalian and avian lineages.

The human X and Y chromosomes, like those
of other animals, are thought to have evolved
from an ordinary pair of autosomes (J). The
pseudoautosomal regions at the termini of the X
and Y chromosomes still recombine during
male meiosis, ensuring X-Y nucleotide se­
quence identity there. Elsewhere on the X and
Y chromosomes, however, X-Y recombination
has been suppressed. These nonrecombining
regions of the X and Y chromosomes have
become highly differentiated during evolution,
and only a few X-Y sequence similarities per-

Howard Hughes Medical Institute, Whitehead Insti­
tute, and Department of Biology, Massachusetts In­
stitute of Technology, 9 Cambridge Center, Cam­
bridge, MA 02142, USA.

*Present address: Department of Human Genetics,
University of Chicago, 924 East 57th Street, Chicago,
IL 60637, USA.
tTo whom correspondence should be addressed. E­
mail: dcpage@wi.mit.edu

sist within them. These modem X-Y gene pairs
are the remaining “fossils” where extensive se­
quence identity between ancestral X and Y
chromosomes once existed. The recent discov­
ery of many X-Y genes has made it possible to
examine the entire group to search for patterns
of human sex chromosome evolution. Thus far,
the human sex chromosomes-the best charac­
terized mammalian sex chromosomes-have
been found to contain 19 X-Y gene pairs (2).

We first compared the locations of all 19
pairs of genes on the human X and Y chromo­
somes (Fig. I). We determined the relative
positions of the X-linked genes through radia­
tion hybrid analysis, in many cases confirming
previously published localizations (3). Map po­
sitions of the Y-linked homologs were obtained
principally from the literature (4-6). On the X
chromosome, most of the X-Y genes map to the
short arm, where they are concentrated toward
the distal end. By contrast, the X-Y genes are

were separated by ultracentrifugation, followed by
two-phase partitioning to enrich for plasma mem­
branes, as described previously [T. W. Short, P. Rey­
mond, W. R. Briggs, Plant Physio/. 101, 647 (1993)].

33. Antibodies against NPH1 were previously described
(7). Rabbit polyclonal antisera were raised (22)
against a COOH-terminal NPH3 fusion protein [CBD­
NPH3C2 (see Fig. 3A)]. CBD-NPH3 protein was ex­
pressed from pET34-Ek/LIC in Escherichia coli and
purified according to manufacturer’s instructions
(Novagen, Madison, WI).

34. NPH1-NPH3 interaction was examined in yeast with
the Matchmaker Gal4 II System (Clontech, Palo Alto,
CA). Expression of fusion peptides was verified by
immunoblot analysis (9, 22) with monoclonal anti­
bodies raised against the Gal4 DNA binding domain
(GBD) and Gal4 activation domain (GAD) (Clontech).

35. J. H. Miller, Experiments in Molecular Genetics (Cold
Spring Harbor Laboratory, Plainview, NY, 1972).

36. We thank R. Harper for data in Fig. 1; J. M. Christie
and W.R. Briggs for GBD-NPH1 constructs and NPH1
antisera; D. Randall for production of NPH3 antisera;
the Arabidopsis Biological Resource Center in Colum­
bus, Ohio, for BAC clones and cDNA libraries; and
members of our laboratory for helpful comments on
the manuscript. This work was funded by USDA
National Research Initiative grant 96-35304-3709,
NSF grant MCB-9723124, and University of Missouri
Research Board grant RB96-055.

3 June 1999; accepted 17 September 1999

found as singletons or small clusters throughout
the euchromatic portion of the Y chromosome.
In general, the map order of the X-linked genes
corresponds poorly to that of the Y-linked ho­
mologs. Local exceptions to this rule are pro­
vided by three small gene clusters that are
present on both X and Y chromosomes (Fig. I).

We next measured, for each of the 19 X-Y
gene pairs, synonymous nucleotide divergence
between the X-linked and Y-linked coding re­
gions (7). Because synonymous substitutions
do not alter the encoded protein, they are gen­
erally assumed to be nearly neutral with respect
to selection. The statistic Ks (the estimated
mean number of synonymous substitutions per
synonymous site) is often used to gauge evolu­
tionary time ( 8). In the present context, Ks
values provide a measure of the evolutionary
time that has elapsed since the gene pairs start­
ed differentiating into distinct X and Y forms.
The calculated Ks values are given in Table I,
where gene pairs are listed according to map
order on the X chromosome.

We noted that the 19 Ks values appeared to
cluster into approximately four groups (Fig. 2):
0.94 to 1.25 (group I), 0.52 to 0.58 (group 2),
0.23 to 0.36 (group 3), and 0.05 to 0.12 (group
4). Each X-Y gene pair’s Ks value differed
significantly from those of all gene pairs in
other groups (P :5 0.02). The most striking
observation was that, on the X chromosome,
the four Ks-defined groups of genes are ar­
ranged in an orderly sequence (Fig. 2). X-Y
genes are stratified by age along the length of
the X chromosome. By contrast, on the Y chro­
mosome, the Ks-defined groups appear to be
scrambled (compare Table I and Fig. I).

What might account for the orderly stratifi­
cation of X-Y genes by age on the human X
chromosome? We hypothesize that, during evo-

29 OCTOBER 1999 VOL 286 SCIENCE www.sciencemag.org

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:10:40 UTC�������������

All use subject to https://about.jstor.org/terms

lution, differentiation of the X from the Y chro­
mosome was initiated one region, or stratum, at
a time. Regions were recruited in the order of
their physical position, with stratum I ( contain­
ing the genes of group I) having been the first
to embark on X-Y differentiation, and stratum 4
having been the most recent. Genes in the same
stratum began differentiating into X and Y ho­
mologs at about the same time, accounting for
their similar Ks values.

X-Y differentiation would have occurred
only after X-Y recombination ceased (9). Our
findings suggest that during evolution, X-Y re­
combination was suppressed regionally, begin­
ning with stratum I and subsequently expanding
in discrete steps to include strata 2, 3, and 4.
Chromosomal inversions, which are known to be
capable of suppressing recombination across
broad regions in mammals ( J 0), would appear to
be the most likely mechanism. These inversions
must have occurred on the evolving Y chromo­
some, where the strata have been scrambled, but
not on the X chromosome, where the order of
strata apparently has been preserved (Figs. I and
2). [Had the strata on the human X chromosome
been extensively shuffled during evolution-as
may have occurred on the mouse X chromosome
after divergence of the human and murine lin­
eages (J J)-we would have observed no corre­
lation between the age ofX-Y gene pairs and the
map positions of their X-chromosomal mem­
bers.] In the modern human sex chromosomes,
the proximal boundary of the pseudoautosomal
region is spanned by a gene that is intact on the
X chromosome, but grossly interrupted on the Y
chromosome (12), consistent with disruption of

REPORTS

an ancient pseudoautosomal region by a Y-chro­
mosomal inversion. We speculate that this par­
ticular event was the most recent in a series of
inversions, each of which enabled X-Y differen­
tiation to begin in one stratum.

This model of staged, region-by-region ini­
tiation of X-Y differentiation also accounts for
two global features of the X chromosome’s
gene content: (i) the concentration in strata 3
and 4 of genes with detectable Y homologs
(Fig. I) and (ii) the concentration on the short
arm (strata 2, 3, and 4) of genes that escape X
inactivation, some with and some without Y
homologs (13). Evolutionary theory predicts
that once X-Y recombination ceased within a
stratum, the genes on the affected portion of the
Y chromosome began to decay, with most of
the Y-linked genes ultimately being obliterated
(J). As an adaptive response, homologous
genes on the X chromosome were up-regulated,
and subsequently became subject to X inactiva­
tion, processes thought to have spread during
evolution on a gene-by-gene or cluster-by-clus­
ter basis (14). If decay of Y-linked genes and
adaptation of X-linked homologs were gradual
evolutionary processes, then one would expect
the youngest X strata to exhibit the highest
densities of (i) genes with detectable Y ho­
mologs and (ii) genes that escape inactivation.
Both predictions are met (Fig. I) (13).

A comparison of the youngest (group 4)
gene pairs with the older (groups I through 3)
gene pairs illustrates certain temporal features of
X-Y differentiation. We measured both synon­
ymous and nonsynonymous substitutions for
each gene pair (Table I). Nonsynonymous sub-

Table 1. Sequence divergence between homologous X- and Y-linked genes.

DNA Protein Sequence
Gene pair Ks KA Ks/KA divergence divergence compared

(%) (%) (nucleotides)

Group 4
GYG2/GYG2P* 0.11 0.06 1.8 7 12 525
ARSDIARSDP* 0.09 0.07 1.3 7 13 846
ARSEI ARSEP* 0.05 0.04 1.2 4 9 615
PRKXIY 0.07 0.03 2.3 5 8 1020
STSISTSP* 0.12 0.10 1.2 11 18 852
KAL1/KALP* 0.07 0.06 1.2 6 12 1302
AMELXIY 0.07 0.07 1.0 7 12 576

Group 3
TB4XIY 0.29 0.04 7.3 7 7 135
EIF7AX/Y 0.32 0.01 32 9 2 432
ZFXIY 0.23 0.04 5.8 7 7 2394
DFFRXIY 0.33 0.05 6.6 11 9 7671
DBXIY 0.36 0.04 9.0 12 9 1932
CASK/CASKP* 0.24 0.22 1.1 15 32 156
UTXIY 0.26 0.08 3.3 12 15 4068

Group 2
UBE7XIY 0.58 0.07 8.3 16 13 693
SMCXIY 0.52 0.08 6.5 17 15 4623

Group 7
RPS4XIY 0.97 0.05 19 18 18 792
RBMXIY 0.94 0.25 3.8 29 38 1188
SOX3/SRY 1.25 0.19 6.6 28 29 264

*Y copy is pseudogene. DNA ‘and protein divergence refer to uncorrected nucleotide (coding region) and amino acid
divergence (nonidentity).

stitutions alter the encoded protein and are con­
strained by selection. Thus, their frequency (KA,
the estimated mean number of nonsynonymous
substitutions per nonsynonymous site) is a func­
tion of both evolutionary time and selective
constraints on the encoded proteins. The degree

X

2

1

GYG2]
ARSD a
ARSE
PRKX
STS ]b
KAL1
AMELX
TB4X
E/F1AX
ZFX

DFFRX] DBX
CASK c
UTX
UBE1X

SMCX

RPS4X

SRY
RPS4Y

ZFY

y

PRKY
AMELY

[ARSEP
a ARSDP

GYG2P

[DFFRY
OBY

c CASKP
UTY

TB4Y

b[ KALP STSP·

SMCY

EIF1AY

RBMY

Fig. 1. Map of homologous
genes in nonrecombining re­
gions of human X and Y chro-

RBMX mosomes. Pseudoautosomal re-
SOX3 gions of X and Y are black; het­

erochromatic region of Y is gray.
Radiation hybrid analysis (3) was
used to map genes on the X
chromosome, which is drawn on
a centiRay scale. Ks-defined stra-

ta on the X chromosome are indicated. The
boundary between strata 2 and 1 is somewhere
between SMCX and RPS4X; here, it is arbitrarily
shown at the centromere (white oval). Genes and
pseudogenes on the Y chromosome were ordered
previously by analysis of naturally occurring de­
letions (4, 5). UBE1 X has a homolog on the squir­
rel monkey Y chromosome but not on the human
Y chromosome (29). Brackets denote three small
gene clusters {labeled a, b, c) that are present on
both X and Y chromosomes.

www.sciencemag.org SCIENCE VOL 286 29 OCTOBER 1999 965

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:10:40 UTC�������������

All use subject to https://about.jstor.org/terms

966

1.25

Group 1

1.00

0 .75

Ks
Group2

0.50

0.25
Group3

p cen q
X ch romosome map position

Fig. 2. Plot of Ks (Table 1) versus X-chromosome
map position (Fig. 1) for 19 X-Y gene pairs.

of constraint can be reflected in the ratio Ksf KA;
values greater than one indicate the presence of
constraints on both homologs, and values in the
vicinity of one are consistent with lack of con­
straint on at least one homolog (8, 15). In groups
1 through 3, 10 of 11 gene pairs exhibit KsfKA
ratios of 3 or higher (Table 1), suggesting that
natural selection has preserved the Y copies of
these genes . Without such selection, these X-Y
homologies (especially those in groups 1 and 2)
would no longer be visible. By contrast, the
seven gene pairs in group 4 show KsfKA ratios
of 1 to 2, and in five of these pairs, the Y copy
is known to be a pseudogene . Among the group
4 pairs, X-Y homology is readily apparent even
in the absence of selective constraint, because
there has been little time for erosion of sequence
similarity. Thus, the Y -chromosomal genes of
the older groups, and especially those of groups

300

~ .. .,
250 >-

15
“‘ C:

200
~
g
.,

150
E . ., .,
g 100 .,
e .,
.2: so ‘O
:>;-
X

0 .08 0 .99 Ks

/
Fig. 3. Plot of X-Y divergence time (age) versus
average Ks value for X-Y gene pairs (weight­
averaged) in each stratum. The X chromosome
schematic is adapted from Fig. 1. Maximum and
minimum age estimates for strata 2, 3, and 4 are
bracketed; these are not statistical confidence
intervals. Theory predicts an approximately linear
relationship between age and Ks value (8); the
shaded area is calibrated with respect to stratum
2, whose age is 130 to 170 million years (21) and
whose average Ks value is 0.53. By extrapolation,
the age of stratum 1 is estimated between 240
and 320 million years.

REPORTS

1 and 2, are survivors of an early winnowing
process that is still ongoing in group 4.

To determine the age of the Ks-defined stra­
ta, we used two methods . First, we considered
published information on homologs of represen­
tative genes in diverse mammals. The maximum
age of stratum 4, for example, was suggested by
the prior observation that homologs of STS and
KALI are pseudoautosomal or autosomal in pro­
simians (16-18). Assuming that suppression of
X-Y recombination is an irreversible evolution­
ary step (14) , this implies that X-Y differentia­
tion in stratum 4 began less than 50 million
years ago (Ma), when the simian and prosimian
lineages diverged (J 9). Minimum ages of the
strata could also be inferred. For example, STS
and KALI have been shown to have X- and
Y-specific homologs in both New and Old
World monkeys (16, 17), suggesting that X-Y
differentiation in stratum 4 began at least 30 Ma,
when the New and Old World monkey lineages
diverged (19, 20). Using similar logic, we in­
ferred the ages of stratum 3 (80 to 130 million
years), stratum 2 (130 to 170 million years), and
stratum 1 (130 to 350 million years) from prior
data on gene homologs in more-distantly related
species, including nonprimate mammals, mar­
supials, monotremes , and birds (21) .

These cross-species comparisons yielded
reasonably precise estimates of age for strata 2,
3, and 4-the younger strata-but only crude
estimates of age for stratum 1. Because this
oldest stratum might contain information about
the origins of mammalian sex chromosomes, its
age is of great interest. Here, we used a second

i,,’
!-;;;-
Ii~

f&
eJ

ti~ c:-
~”” -~ l§ 1:!”ii,-

~~

dating method, based on Ks values for X-Y gene
pairs. Theory predicts that among human X-Y
gene pairs, Ks values should be roughly propor­
tional to age ( 8). This expectation is met by the
X-Y gene pairs of strata 2, 3, and 4 (Fig . 3). By
extrapolation, we estimated that X-Y differenti­
ation began 240 to 320 Ma in stratum 1 (Fig . 3).
These findings suggest that X-Y divergence be­
gan shortly after the mammalian lineage arose,
having diverged from the lineage of birds (with
Z-W sex chromosomes) between 300 and 350
Ma (19). [Because the sex chromosomes of
birds appear to be completely unrelated to the
mammalian sex chromosomes, it is thought that
they arose independently, from a different auto­
somal pair (22).] Interestingly, our Ks findings
indicate that SOX3 and SRY (the primary sex­
determining gene) are among the oldest known
X-Y gene pairs in humans (Table 1). This find­
ing strengthens an hypothesis, by Foster and
Graves , which states that an ordinary autosomal
pair became sex chromosomes when mutations
fashioned one allele of SOX3, originally an au­
tosomal gene, into the male-determining factor
SRY (23). Indeed, formal cluster analysis of the
Ks values we report suggests that the X-Y genes
of group 1 might actually comprise two distinct
strata, with SRY!SOX3 perhaps being older than
the two other X-Y gene pairs of group 1
(RPS4XIY and RBMXIY) (24). Although the dif­
ference in Ks values between SRY!SOX3 and the
two other X-Y gene pairs is not statistically
significant, the evidence is suggestive .

If future studies establish that the group 1
genes are divisible into two strata, these results

~
[
~ .s ;–

l/ 4
~ …. c:- .§ tf , .J2 3 3 q~ 1:!;;;- G-=—?i°#

fR ?i = ~~ ::..$ 0~ 2 g .f 2 ·S’O 2 ·S 2 ~c:- bi ii’ .::… ::? .::… fil ~ § ~~ ,!::’ ~-;;; (J …. ~ §~
iff~ 0 ‘- i1f ~= if Clj

Autosome X y X y X y X y XY

! ! ! !
humans

Autosome XYin XYin XY in non-simian
in birds monotremes marsupials placental mammals

Fig. 4. A proposed sequence of evolutionary events that generated four strata on the human X
chromosome. Four inversions on the Y chromosome are postulated . Each inversion reduced the size
of the pseudoautosomal ( X-Y recombining) region (black; for simplicity, only one pseudoautosomal
region is shown for each chromosome) and enlarged the portions of the X (yellow) and Y (blue)
chromosomes that did not recombine during male meiosis. Ongoing decay and loss of Y genes
offset these periodic expansions of the nonrecombining region of the Y chromosome. Points of
divergence from the sex chromosomes of other mammals are indicated. This model does not
preclude the occurrence of (i) additional inversions or other rearrangements within the nonrecom ­
bining portion of the evolving Y chromosome or (ii) similar rearrangements on the evolving X
chromosome, so long as they do not disturb the fundamental order among the four strata.

29 OCTOBER 1999 VOL 286 SCIENCE www .sciencemag.org

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:10:40 UTC�������������

All use subject to https://about.jstor.org/terms

would also help date the emergence of X inac­
tivation during mammalian sex chromosome
evolution. XIST, an X-specific gene which plays
a pivotal role in X inactivation (25), is located
near RPS4X and therefore would be in the
younger of the two strata~not in the stratum
where the nascent X and Y chromosomes first
differentiated. This would controvert the hypoth­
esis of Chandra, who speculated that X inactiva­
tion emerged contemporaneously with the chro­
mosomal sex-determining mechanism (26).

Consistent with our evolutionary map,
Graves and colleagues have postulated that the
long arm and proximal short arm of the human
X chromosome are at least 170 million years
old (27, 28). They have referred to this portion
of the X as the “XCR” (X conserved region).
Graves’s XCR corresponds approximately to
our strata 1 and 2. They have also postulated
that the distal short arm of the human X chro­
mosome is younger. This “XAR” (X added
region) was attributed to translocation of an
autosome to the pseudoautosomal region of
both X and Y after divergence of placental
mammals from marsupials (27, 28). Our strata
3 and 4 are found within Graves’s XAR.

In conclusion, we postulate that the evolution
of human sex chromosomes was punctuated by
at least four events, plausibly a series of inver­
sions on the Y chromosome (Fig. 4). Each event
suppressed X-Y recombination in one stratum
and enabled X-Y differentiation to proceed
there. The first of these events, which created
stratum 1, was roughly contemporaneous with
the birth of the mammalian sex chromosomes
and the emergence of SRY as the primary sex
determinant. This occurred about 240 to 320 Ma,
shortly after the mammalian and avian lineages
diverged. The pseudoautosomal region was ex­
panded by translocation of autosomal material
between the second and third events (which cre­
ated strata 2 and 3, respectively). The fourth
event occurred relatively recently, during pri­
mate evolution, creating stratum 4, where X-Y
differentiation is still in its earliest stages.

References and Notes
1. J. J. Bull, Evolution of Sex Determining Mechanisms (Ben­

jamin Cummings, Menlo Park, CA, 1983); J. A. Graves,
Annu. Rev. Genet. 30, 233 (1996); B. Charlesworth,
Curr. Biol. 6, 149 (1996); W.R. Rice, Bioscience 46, 331
(1996).

2. The 19 X-Y gene pairs studied include the following:
GYG2/GYG2P [J. Mu, A. V. Skurat, P. j. Roach,}. Biol.
Chem. 272, 27589 (1997); (6)], ARSDIARSDP, ARSE/
ARSEP [G. Meroni et al., Hum. Mo/. Genet. 5, 423
(1996)], PRKX!Y[A. Klink et al., Hum. Mal. Genet. 4,869
(1995); K. Schiebel et al., Hum. Mal. Genet. 6, 1985
(1997)], STSISTSP (16), KAL1/KALP [B. Franco et al.,
Nature 353, 529 (1991); R. Legouis et al., Cell 67,423
(1991); (77)], AMELXIY [Y. Nakahori, 0. Takenaka, Y.
Nakagome, Genomics 9,264 (1991)], TB4XIY [H. Gonda
et al.,}. /mmunol. 139, 3840 (1987); (5)], ZFXIY [D. C.
Page et al., Cell 51, 1091 (1987); A. Schneider-Gadicke,
P. Beer-Romero, L G. Brown, R. Nussbaum, D. C. Page,
Cell 57, 1247 (1989)], EIF1AXIY [T. E. Dever et al.,}. Biol.
Chem. 269, 3212 (1994); (5)], DFFRXIY [M. H. Jones et
al., Hum. Mal. Genet. 5, 1695 (1996); (5)], DBXIY (5),
CASK/CASKP [A. R. Cohen et al.,}. Cell Biol. 142, 129
(1998); (6)], UTX/Y(5),SMCXIY[J. Wu eta/., Hum. Mal.

REPORTS

Genet. 3, 153 (1994); A. I. Agulnik et al., Hum. Mal.
Genet. 3, 879 (1994)]. RPS4XIY [E. M. Fisher et al., Cell
63, 1205 (1990)]. RBMXIY [M. Soulard et al., Nucleic
Acids Res. 21, 4210 (1993); K. Ma et al., Cell 75, 1287
(1993); M. L. Delbridge, P.A. Lingenfelter, C. M. Disteche,
J. A. Graves, Nature Genet. 22, 223 (1999); S. Mazeyrat,
N. Saut, M. G. Mattei, M. J. Mitchell, Nature Genet. 22,
224 (1999)], SOX3/SRY [M. Stevanovic, R. Lovell-Badge,
J. Collignon, P. N. Goodfellow, Hum. Mo/. Genet. 2,
2013 (1993); A. H. Sinclair et al., Nature 346, 240
(1990)]. One interspecies pair was also studied: human
UBE1 X [P. M. Handley, M. Mueckler, N. R. Siegel, A.
Ciechanover, A. L. Schwartz, Proc. Natl. Acad. Sci. U.S.A.
88, 258 (1991)] and squirrel monkey UBE1Y (29). In
humans, UBE1Y was deleted from the Y chromosome
(29). We used squirrel monkey UBE1 Y as a substitute.

3. Using polymerase chain reaction (PCR), we tested DNAs
from the 93 hybrid cell lines of the GeneBridge 4 panel
(Research Genetics) [G. Gyapay et al., Hum. Mo/. Genet.
5, 339 (1996)] for the presence of each of the X-linked
genes. PCR conditions and primer sequences have been
deposited at GenBank, where accession numbers are as
follows: GYG2, G49430; ARSD, G42687; ARSE, G42688;
PRKX, G42689; STS, G42690; KAL1, G42691; AMELX,
G42692; TB4X, G34979; E/F1AX, G34989; ZFX, G42693;
DFFRX, G34982; DBX, G34988; CASK, G49441; UTX,
G34976; UBE1X, G42694; SMCX, G42695; RPS4X,
AF041428; RBMX, G42696; and SOX3, G42697. Analy­
sis of the results positioned the genes with respect to
the radiation hybrid map of the X chromosome con­
structed at the Whitehead/MIT Center for Genome
Research [T. J. Hudson et al., Science 270, 1945 (1995);
www-genome.wi.mit.edu/cgi-bin/contig/phys_map].

4. D. Vollrath et al., Science 258, 52 (1992).
5. B. T. Lahn and D. C. Page, Science 278, 675 (1997).
6. C. Sun et al., Nature Genet., in press.
7. Homologous X and Y DNA sequences were aligned by

means of MegAlign software (DNASTAR, Madison, WI).
For each X-Y gene pair, estimates of the mean numbers
of synonymous substitutions per synonymous site (Ks),
and of nonsynonymous substitutions per nonsynony­
mous site (KA)-all corrected for multiple changes-were
calculated using published algorithms (8) as implemented
in GCG software (Genetics Computer Group, Madison,
WI). Insertions and deletions were ignored in these cal­
culations. In the case of SOX3 and SRY, sequence similar­
ity is limited to, and our analysis was restricted to, the
HMG box domain. Our analyses of other X-Y gene pairs
employed all available coding sequences. Only a partial
UBE1Y (squirrel monkey) coding sequence was available
for comparison with its human X homolog. Sequences for
all pseudogenes were extracted from genomic sequences:
GYG2P, ARSDP, and ARSEP from BAC (bacterial artificial
chromosome) clone 203M13 (GenBank AC002992); STSP
from BAC clone NH0494J04 (GenBank AC006382); KALP
from BAC clone NH0292P09 (GenBank AC006370);
CASKP from BAC clone 47511 (GenBank AC004474). Se­
quences for all other genes were obtained from published
cDNAs, whose GenBank accession numbers are as fol­
lows: GYG2, U94362; ARSD, X83572; ARSE, X83573;
PRKX, X85545; PRKY, Y15801; STS, M16505; KAL1,
M97252; AMELX, M86932; AMELY, M86933; TB4X,
Ml 7733; TB4Y, AF000989; ZFX, X59739; ZFY, M30607;
EIF1AX, L 18960; E/F1AY, AF000987; DFFRX, X98296;
DFFRY, AF000986; DBX, AF000982; DBY, AF000984;
CASK, AF032119; UTX, AF000992; UTY, AF000994;
UBE1X, M58028; UBE1Y, AJ003105; SMCX, L25270;
SMCY, U52191; RPS4X, M58458; RPS4Y, M58459; RBMX,
Z23064; RBMY, X76059; SOX

Science
~MAS

D
JSTOR

Natural Selection and the Origin of jingwei, a Chimeric Processed Functional Gene in
Drosophila

Author(s): Manyuan Long and Charles H. Langley

Source: Science , Apr. 2, 1993, New Series, Vol. 260, No. 5104 (Apr. 2, 1993), pp. 91-95

Published by: American Association for the Advancement of Science

Stable URL: https://www.jstor.org/stable/2881129

JSTOR is a not-for-profit service that helps scholars, researchers, and students discover, use, and build upon a wide
range of content in a trusted digital archive. We use information technology and tools to increase productivity and
facilitate new forms of scholarship. For more information about JSTOR, please contact support@jstor.org.

Your use of the JSTOR archive indicates your acceptance of the Terms & Conditions of Use, available at

Terms and Conditions of Use

American Association for the Advancement of Science is collaborating with JSTOR to digitize,
preserve and extend access to Science

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:13:06 UTC�������������

All use subject to https://about.jstor.org/terms

111:~~I HlllllllH’lllf-‘lllillllllllllll®I-Hll11!11&111elllll!llill11flll!llll!!IIH Jll■l■lallll REPORTS

Phosphoinositides could serve as specific
membrane targets that bind proteins re­
quired for the formation of transport vesi­
cles, such as the polypeptides of the adaptor
complex that link clathrin to the cytoplas­
mic tails of certain transmembrane receptor
proteins {for example, the mannose-6-phos­
phate receptor) (28). Vesicles from rat adi­
pocytes that contain the glucose transporter
also contain PI 4-kinase, which may regulate
the transport (fusion) of these vesicles with
the plasma membrane in response to insulin
(29). In addition to its potential role in
signaling cell proliferation, PI 3-kinase asso­
ciated with receptor protein tyrosine kinases
at the plasma membrane may also take part
in the endocytosis and down-regulation {ly­
sosomal degradation) of these receptors. In
this way, the duration and the magnitude of
the growth signal might be modulated. The
association of Vps34p with the membrane
appears to be mediated by the product of
another VPS gene, VPS15. The VPS15 gene
encodes a membrane-associated protein ki­
nase (Vpsl5) (30, 31) that can be chem­
ically cross-linked to Vps34p (21). This
raises the possibility that the Vpsl5 and
Vps34 proteins may function together as
components of a signal transduction com­
plex that regulates intracellular protein
sorting decisions.

REFERENCES AND NOTES

1. L. C. Cantley et al., Cell 64, 281 (1991).
2. J. M. Backer et al., EMBO J. 11, 3469 (1992).
3. S. Soltoff, S. Rabin, L. Cantley, D. Kaplan, J. Biol.

Chem. 267, 17472 (1992).
4. M. J. Berridge and R. F. Irvine, Nature 341,197

(1989).
5. M. Whitman et al., ibid. 332,644 (1988).
6. D. L. Lips et al., J. Biol. Chem. 264, 8759 (1989).
7. L.A. Serunian et al., ibid., p. 17809.
8. C. Carpenter et al., ibid. 265, 19704 (1990).
9. F. Shibasaki, Y. Homma, T. Takenawa, ibid. 266,

8108 (1991).
10. M. Otsu et al., Cell 65, 91 (1991 ).
11. I. D. Hiles et al., ibid. 70,419 (1992).
12. C. A. Koch, D. Anderson, M. F. Moran, C. Ellis, T.

Pawson, Science 252, 668 (1991).
13. P. Hu et al., Mo/. Cell. Biol. 12, 981 (1992).
14. C. J. McGlade et al., ibid., p. 991.
15. P. K. Herman and S. D. Emr, ibid. 10, 6742 (1990).
16. S. Kornfeld and I. Mellman, Annu. Rev. Cell Biol.

5, 483 (1989).
17. D. J. Klionsky, P. K. Herman, S. D. Emr, Microbial.

Rev. 54, 266 (1990).
18. J. E. Rothman and L. Orci, Nature 355, 409

(1992).
19. T. Stevens et al., Cell 30, 439 (1982).
20. T. R. Graham and S. D. Emr, J. Cell Biol. 114, 207

(1991).
21. J. H. Stack, P. K. Herman, P. V. Schu, S. D. Emr,

EMBO J., in press.
22. K. R. Auger, C. L. Carpenter, L. C. Cantley, L.

Varticovski, J. Biol. Chem. 264, 20181 (1989).
23. G. Endemann, S. N. Dunn, L. C. Cantley, Bio­

chemistry26, 6845 (1987).
24. Genetics Computer Group sequence analysis

package; University of Wisconsin.
25. S. K. Hanks eta/., Science 241, 42 (1988).
26. D. R. Knighton et al., ibid. 253, 407 (1991).
27. M. P. Sheetz and S. J. Singer, Proc. Natl. Acad.

Sci. U.S.A. 71, 4457 (1974).
28. B. M. Pearse and M. S. Robinson, Annu. Rev. Cell

Biol. 6, 151 (1990).

29. R. L. Del Vecchio and P. F. Pilch, J. Biol. Chem.
266, 13278 (1991).

30. P. K. Herman et al., Cell 64, 425 (1991).
31. P. K. Herman, J. H. Stack, S. D. Emr, EMBOJ.10,

4049 (1991).
32. T. Kunkel, Proc. Natl. Acad. Sci. U.S.A. 82, 5463

(1985).
33. F. Sherman, G. R. Fink, L. W. Lawrence, Methods in

Yeast Genetics: A Laboratory Manual (Cold Spring
Harbor Laboratory, Cold Spring Harbor, NY, 1979).

34. M. Whitman et al., Nature 315,239 (1985).
35. J. P. Walsh, K. K. Caldwell, P. W. Majerus, Proc.

Natl. Acad. Sci. U.S.A. 88, 9184 (1991).
36. L. J. Wickerham, J. Bacterial. 52, 293 (1946).
37. K. R. Auger et al., Cell 57, 167 (1989).
38. We thank members of the Emr lab for discussions

and G. Huyer” for his help in optimizing the Pl
3-kinase assay. Supported by the Howard Hughes
Medical Institute, a grant from the NSF (to S.D.E.),
and by a fellowship from the Deutsche Forschungs­
gemeinschaft (to P.V.S.). S.D.E. is an investigator of
the Howard Hughes Medical Institute.

8 October 1992; accepted 13 January 1993

Natural Selection and the Origin of jingwei, a
Chimeric Processed Functional Gene in Drosophila

Manyuan Long* and Charles H. Langley

The origin of new genes includes both the initial molecular events and subsequent pop­
ulation dynamics. A processed Drosophila alcohol dehydrogenase (Adh) gene, previously
thought to be a pseudogene, provided an opportunity to examine the two phases of the
origin of a new gene. The sequence of the processed Adh messenger RNA became part
of a new functional gene by capturing several upstream exons and introns of an unrelated
gene. This novel chimeric gene, jingwei, differs from its parent Adh gene in both its pattern
of expression and rate of molecular evolution. Natural selection participated in the origin
and subsequent evolution of this gene.

How genes with novel functions evolve
remains a fundamental and fascinating ques­
tion. Gene duplications (1), exon shuffling,
and processed genes (2) have been suggested
as important sources of novel genes, but
little is known about the evolutionary mech­
anisms or the participation of natural selec­
tion in their early history. Our analysis of
the structure, expression, and evolution of a
putative processed Adh pseudogene in Dro­
sophila provided an opportunity to examine
both the early molecular events and the
evolutionary processes that created it. The
jingwei (igw) gene, as we named it (3), is a
locus located on chromosome 3 in the Dro­
sophila sibling species D. teissieri and D.
yakuba. A part of jgw was initially observed
to hybridize to the Adh probe (4). Further
analysis suggested that in a single event, in
the ancestor of the two species the Adh
portion of this potential pseudogene was
retrotransposed from an mRNA of the Adh
locus on chromosome 2 (5). To understand
the molecular population genetics of this
potential pseudogene, we investigated its
within-species DNA sequence variation.
However, our results convince us thatjgw is
not a pseudogene and that the Adh-derived
sequence is a part of a novel gene.

The Adh-derived portion was sequenced
from ten jgw alleles of D. teissieri and of 20
jgw alleles of D. yakuba collected from
natural populations (6). Figure 1 shows the

Section of Genetics, Section of Evolution and Ecology,
and Center for Population Biology, University of Cali­
fornia at Davis, Davis, CA 95616.

•To whom correspondence should be addressed.

SCIENCE • VOL. 260 • 2 APRIL 1993

distribution of nucleotide polymorphisms
within species and the variation between
species in a 765-bp segment that corre­
sponds to the protein coding region of the
Adh gene. Only one possible ancestral poly­
morphism was apparent {site 782). A sum­
mary of the DNA sequence variation in this
region is presented {Table 1).

Unexpectedly, most polymorphisms were
silent {eight out of ten in D. teissieri and 19
out of 21 in D. yakuba). If jgw were a pseudo­
gene in which mutations had no phenotypic
effect, most changes would be at replacement
sites and there would be a reasonable frequen­
cy of stop codons (7). No new stop codons
were present. Another prediction of the pseu­
dogene hypothesis is that a pseudogene has a
higher overall level of nucleotide variation.
Estimates of DNA sequence polymorphism in
jgw were similar to that found in the Adh gene
{Table 1) and are typical of many functional
genes in Drosophila (=0.005) (8). A compar­
ison of between-species divergence also re­
vealed a significant bias toward silent versus
replacement substitutions, and the degree of
bias was smaller than for the within-species
comparison. The bias was noted in the earlier
study (5) but by itself was not large enough to
convince the authors that this was not a
pseudogene. In addition, insertions and dele­
tions are abundant in the evolution of mam­
malian pseudogenes (7). In jgw, no length
polymorphism or divergence was observed in
the coding region (9). These molecular pop­
ulation genetic observations imply that the
Adh-derived sequence is all or part of jgw, a
new functional gene in D. teissieri and D.
yakuba. This conclusion motivated our search

91

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:13:06 UTC�������������

All use subject to https://about.jstor.org/terms

D.teissieri D.yalcuba Fig. 1. Polymorphic
and divergent differ­
ences in the Adh-de­
rived portion of jgw.
The nucleotide posi­
tions are indicated in
the first column, num­
bered according to that
reported from D. yaku­
ba (5). The second and
thirteenth columns show
the consensus nucleo­
tides. The dots indicate
that the nucleotides are
identical to those of the
consensus. The first row
contains the numbers of
the sequenced alleles
(10 from D. teissieri and
20 from D. yakuba).

I 2 3 4 5 6 7 I 9 10 I 2 3 • 5 6 l I • 10 11 12 ll 14 IS 16 17 IM 1fl 20
219
225
231
232
2•s
2•6
254
297
323
335
365
369
370
372
3H
387
388
389
390 …
08
•so
•s1
09
•60
•61
•63 .,.
00 ••2
•97
sea
525
526
527
535
557
566
569
581
599
608
614
650
680
719
728
70 ,.,
776
779
782
785
112
863
178
884
887
893

••• 917
9U
951
960
961
974
977

G
A
T
C
G
C
G
C
G
C
G
G
C
C
C
T
G
C
Q
G
A
A
C
Q
Q
C
A
A
Q
C
C
T
C
G
C
C
Q
G
C
G
G
G
Q
C
C
G
C
G
T
C
C
A
T
T
Q
Q
C
Q
Q
C
A
C
A
T
Q
C
G

i: ;.

;.

o

o

;-

for an RNA from jgw and a more detailed
analysis of its structure.

Northern (RNA) analysis (10) of both
total and polyadenylate [poly(A) ]-selected
RNAs from adults probed with jgw yielded
only a single band that corresponded in size
to that derived from the Adh gene (11).
Because of the sequence similarity between
the two genes and the abundance of the Adh
mRNA, this result was not surprising. Two
interpretations are possible. (i) The jgw

A

;- ;- ;-
;.

A i:: A C

G
T
A
G
A
A
A
G
G
C
G
A
A
A
A
G
C
T
A
T
G
T
G
T
C
C
G
G
A
C
T
C
G
C
A
C
Q
Q
Q
A
G
C
T
C
C
A
C
G
T
C
C
A
G
T

o

;.

o
C i::

o

o o o

G G

Q T
Q A
T
Q
C
C
C
A
A
A
C
T
Q

f

;-

o

;-

;.

;-

mRNA may be detectably abundant but
obscured by the Adh signal if its size is similar
to that of Adh. (ii) The jgw mRNA may be
so rare (or nonexistent) that it was undetect­
able by the Northern technique.

We used the polymerase chain reaction
(PCR) in conjunction with reverse transcrip­
tion (RT-PCR) (12) to identify RNAs with
greater sensitivity and specificity than is pos­
sible with Northern analysis. Several regions
can be found in which Adh and the Adh-

Table 1. Summary of the variation within and between species. The numbers of segregating (or
polymorphic) sites are classified as replacement (R) and silent (S). Estimates of 8 and ‘II’ are
measures of within-population variation in terms of the parameter 4Nµ per nucleotide, assuming the
neutral theory of molecular evolution (20) where N is the effective population size and µ is the
mutation rate for selectively equivalent nucleotides. The values of divergence between species
(average divergence per nucleotide) are the results of averaging pairwise comparisons between
the alleles from the two species minus average within-species differences, subject to the multiple
substitution correction according to the Jukes and Cantor model (28). The standard deviations of
8 and divergence estimates were calculated as described (29). Adh polymorphism data (19 of all
the samples) of 0. yakuba and D. teissieri Adh sequences are from (5, 19).

Gene

D. teissieri jgw
D. yakuba jgw
D. yakuba Adh

jgw
Adh

92

Segregating sites (n)

R s
Polymorphism within species

2 8
2 19
0 18

Divergence between species
0.041 ± 0.009 0.127 ± 0.027
0.007 ± 0.004 0.059 ± 0.018

8

0.005 ± 0.003
0.008 ± 0.003
0.006 ± 0.003

‘II’

0.005
0.006
0.006

SCIENCE • VOL. 260 • 2 APRIL 1993

f

;.

;.

;.

C

;.

C

A A

A A A A

(i 0
Ci (i

i T

i:
! A 1 ! T

i. A A A A
T T T

f

;.

G 0

i ;.
;.

T

derived portion of jgw differ by two or three
contiguous substitutions. Under optimal PCR
conditions, primers that contain these substi­
tutions at their 3′ ends specifically amplified
jgw RNA. A jgw-specific amplification prod­
uct was obtained from total RNA extracted
from both species. Sequencing these products
confirmed the characteristic substitutions of
the jgw gene in the 5 71-bp amplified fragment
from D. teissieri and the 321-bp amplified
fragment from D. yakuba (13). These results
demonstrated the presence of jgw-derived
RNA in both species.

The single band in the 5′ rapid amplifica­
tion cDNA end (RACE) product and the
identity of the two sequenced clones from this
product indicated thatjgw RNA has a distinct
5′ end in D. teissieri (14). In both species, jgw
appears to have captured more than 180 hp of
additional exon (or exons) in the 5′ direction
(Fig. 2). The length of the 3′ end ofD. teissieri
jgw RNA is similar to that of Adh RNA.
Figure 2 also shows the sequence of these
products. By extending the Adh reading frame
in the 5′ direction, the D. teissieri jgw appears
to have acquired a new start codon and
potentially encodes a hybrid protein with an
additional 77 amino acids added to the Adh­
derived protein. A preliminary analysis of jgw
in D. yakuba indicated that it contains a
similar structure.

To determine the structure of the genomic
region (or regions) containing the exon (or

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:13:06 UTC�������������

All use subject to https://about.jstor.org/terms

Fig. 2. Structure of the
jgw transcription unit
(13). (A) The structure
of jgw deduced from
the comparison of the
sequenced 0. teissieri
RNA and genomic
DNA The four boxes
show the four exons in
jgw , and the solid and
hatched portions of the
boxes are the putative
protein coding regions,

A
+-Captured region-j

ATG

/lngwel

Adh

B D. telsslerl 5′ portion sequence of jgw
AGTATCAAAGT’M’CATTGTATTAGAACAAACATM”CAA’J’TTCGGC.MACATTACAAAAT.v.AAAAAMAATC.AAAAAATTfA

CTGTATTGTATTAACCAA’l”M’GA’ffATAUSiGCGCTTCGTCTTACAACCATTACTCTTCTAAAACGCACTCCTTTATTACTCGTGCCAAAGg
Exon 1

t9a9tctacaaaccctattag:atcctgttcatat99taattcattattaacaat9c9attc9cagTTAAAAGCAGCACCGACMCACGCACT
lntron 1 Exon 2

CTACTACTAAGAAgt.Aagtaactcgt.9ggeaaH9tccaa9tca9aatt9ctut9t9caaacecatGTCCTTTTCCAGGGAATCTGGTACC
lntn>n2

GTTTTCMa.TGATGA~TCATGCM:CATGGGGGGgtaagtactatatgatgagttgtctagaactcca
Exon 3

uactca.agacttaccacct.ttgtcataaacctataaaaggattataaaagatatatatatatatattatttttcagTGCGGTTGTCAGTTG -· ~~~ ……… .
Exon•~

Adh-derived and recruited, respectively. The open boxes are the putative
untranslated portions . The lines among boxes are the three intrans. The
exon structure of the original Adh gene is also given , based on (5) . The
gray portions in four Adh exons are fused into one exon in the jgwgene. (8)
Sequences of the 5′ regions of jgwlrom 0. teissieri and 0 . yakuba. Intrans
are shown in lowercase letters. The codon ATG underlined in the 0.
teissieri sequence is the hypothetical start codon. The 3′ region of 0.
teissieri jgw is the same as that of the published sequence (5). An
apparent polyadenylation signal, AATAAA, is at positions 1123 to 1128 of
0. teissieri.

D. yakube 5′ portion sequence of Jgw
1 …. .. 7’1’C1″‘J’TCGACCA~CTTTGTTACTCGTGCAAAAGgtgaggctaaaaatccctattagattctgttc

Exon 1 lntron 1

atat99taattcactattucaatqecattcgcagTTAAAAGCAGCACAGCCGCACTCTACTACTCTCTACTACMCAAgtaa9taa
Exon 2

ctc9t999caua9tccaa9tca99att,ctu.9t9caaacccatqtccttttacaqGGM.TCTGGT’J’CCGTTTAGGTCATGATTCTTGT
0-.2 Exon3

GTGTTCTCAAATMGGCCATCATCCCCCATGGGGGCg”ta9t-agtcta9uctccaaucttaccacattt9atataa9cgattcaaat9aga

caaatacataacutcca99a9atctgctee99gt.tctcttctcttcatatc9atcautuauttgt.ctca9t9ca9ttgtca9ttgcag
lnln>n3

ttca9C9agUattgcatctctttaaattUcttatttgatcaaatc9aca~ . … …. .
Exon•~

exons) from which this new mRNA is de­
rived, we carried out PCR-based genomic
sequencing. Three intrans (and three exons)
were found 5′ to the Adh-derived exon in
both species (Fig. 2). The size and positions of
the three intrans are similar in the two spe­
cies, and the standard intron-splicing signal
(GT-AG) was found in each case except for
the D. teissieri second intron (GT-AT).

Because the structure of jgw is similar in
both species, it seems likely that the insert­
ed exons of Adh recruited the observed
three 5′ exons (and introns) rather quickly.
This interpretation is supported by the ob­
servation that no silent substitution accu­
mulated in the Adh-derived portion of jgw
before speciation, although several replace­
ment substitutions were common to jgw in
both species. This suggests that there was
relatively little time between the insertion
event and the subsequent molecular evolu­
tionary events: the recruitment of the addi­
tional 5′ components into a functional
transcription unit, the substitution of sev­
eral amino acids, and the speciation.

In spite of similarities of structure and early
common ancestry, jgw appears to have
evolved distinct patterns of expression in the
two species. Figure 3 shows the pattern of jgw
expression at different development stages de­
termined by quantitative PCR technique
( I 5). Quantitative PCR provides at least a
relative comparison of abundance and con­
firmed that the amount of jgw RNA was
generally smaller than the amount of Adh. In
D. teissieri, jgw showed sex-specific expression.
The expected band was abundant in amplifi­
cations of RNA from adult males, whereas no
band was present in reactions from RNA of
other stages or adult females. However, in D.
yakuba reactions with RNA from all stages
showed less expression than was seen with D.
teissieri RNA. Larval stages L1 and L2 and

adult male RNA showed more expression
than L3 and adult female RNA. Just as the
protein sequences evolved rapidly, the se­
quences controlling the regulation of jgw have
also evolved rapidly. Further analysis will
determine if the species-specific patterns of
expression are the result of the divergence of
sequences in and around the jgw transcription
unit or whether these expression patterns
reflect divergence at more distant, trans-act­
ing regulatory loci.

Whereas a distinct mRNA was expressed
fromjgw in each species, the conclusion that
jgw is a functional gene also rests on the
evolutionary interpretation of the polymor­
phism observed within each species. The
extension of this evolutionary analysis in
order to incorporate between-species diver­
gence reveals that the early history of jgw
was dominated by positive natural selection
instead of the neutral mutation and genetic
drift that are expected to characterize the
dynamics of pseudogenes (16).

To explore the divergence of jgw in the
period immediately after the Adh retro-

Fig. 3. Developmental pattern of jgwexpression .
The same method as described ( 12) was used
to amplify jgw RNA from total RNA extracted
from the different development stages. The
cDNA product from 300 ng of RNA was added
as a template for amplification because for the
range from 100 ng to 1000 ng of total RNA, the
yield of jgw RNA PCR products was observed to
be linearly related to the template RNA concen­
tration ( 11). PCR cycles (30) were conducted at
94°C (1 min) for denaturation, 60°C (2 min) for
annealing, and 72°C (3 min) for polymerization.
The products of the PCR reactions were electro­
phoresed on a 1 % agarose gel and stained with
ethidium bromide . More Adh RNA was amplified
from the same cDNA preparations, which sug­

transposition but before the split of D.
teissieri and D. yakuba, we compared the
Adh-derived portions of all jgw alleles with
the available Adh alleles in D. teissieri and
D. yakuba and with Adh alleles from various
members of the melanogaster subgroup (Fig.
4). The fixed nucleotide substitutions that
distinguished the Adh-derived portion of all
jgw alleles from Adh in all alleles of D.
teissieri and D. yakuba were assumed to have
been fixed in the ancestral species after the
insertion. Eight such sites exist at positions
339, 576, 600, 791, 831, 861, 918, and
954. Surprisingly, all eight substitutions
were amino acid replacements; no silent
substitutions were found. By incorporating
the available Adh sequences of the other six
members of the melanogaster subgroup, we
found that all eight substitutions were
unique mutations that were fixed in the jgw
lineage. Assuming the replacement rate of
substitution of 0.5 x 10- 9 per nucleotide
per year [calculated from melanogaster sub­
group Adh data (I 7) I, these eight substitu­
tions (out of 572 possible replacement sites)

D. teissieri
Adult

L1 L2 L3 Female Male

571

u

321

gests that Adh is expressed in larger amounts than jgw ( 11). Molecular size markers are on the right
in base pairs .

SCIENCE • VOL. 260 • 2 APRIL 1993 93

This content downloaded from
�������������71.57.134.161 on Wed, 29 Dec 2021 11:13:06 UTC�������������

All use subject to https://about.jstor.org/terms

Fig. 4. Phylogenetic
analysis of the evolution
of the Adh-derived por­
tion of jingwei. The topol­
ogy reflects the as­
sumed evolutionary rela­
tion between jingwei and
Adh (5). The small (stip­
pled) trees at the tips in­
dicate that the detailed
phylogenies of the al­
leles within species are

outgroup
species

teissieri
jingwei
yakuba

unknown. Numbers presented as fractions are the numbers of (fixed or polymorphic) amino acid
replacement differences over the numbers of (fixed or polymorphic) silent differences. The Adh outgroup
includes the sequences of 0. melanogaster, 0. simulans, 0. sechellia, 0. mauritiana, 0. erecta, and 0.
arena. These outgroup sequences were used to confirm that all functional Adh genes have the same
DNA sequence for the eight codons that are distinct to jingwei (in both 0. teissieri and 0. yakuba) with
one exception. At the eighth codon (nucleotide position 954) of 0. arena Adh, there is a substitution that
conferred an amino acid replacement. The italic numbers are the average numbers of silent substitutions
and the numbers of silent polymorphisms under the HKA test, which also incorporates the standard
assumptions of the neutral theory.

yielded an estimate of about 30 million
years for the time of the last common
ancestor of jgw and Adh. This conflicts with
the age of the melanogaster subgroup, which
is estimated to be 17 to 20 million years
(18). The lack of silent substitutions im­
plies that the insertion of jgw occurred close
to the· time of the speciation of D. teissieri
and D. yakuba. The excess of replacement
substitutions in the jgw line of descent
between the time of its insertion and the
speciation of D. teissieri and D. yakuba is
consistent with the model thatjgw respond­
ed to positive natural selection and evolved
a new function.

This analysis of the early history of jgw
leads us to examine its dynamics after the
divergence of the two species. One method to
explore the role of natural selection at a single
locus is to compare the replacement and silent
substitutions found between species with the
replacement and silent polymorphism found
within species ( I 9). The variation in the
Adh-derived portion of jgw within and be­
tween D. teissieri and D. yakuba is summarized
(Fig. 4). A relative excess of (fixed) replace­
ment substitutions over (fixed) silent substitu­
tions between species (21: 16) is apparent
when compared to the proportion of replace­
ment polymorphisms over silent polymor­
phisms within a species (4:27). These results
are inconsistent with the null hypothesis in
the test proposed (19), which incorporates
the assumption of the selective neutrality (x2
= 13.6, P < 0.001). Therefore, adaptive
protein evolution remained an important
force in the evolutionary history of jgw after
the separation of the two species. Further
evidence supporting this view of the diver­
gence of jgw in D. teissieri and D. yakuba
comes from a comparison of the captured 5′
coding regions. Eleven out of the 15 differ­
ences in 153 alignable coding sites cause
amino acid replacements.

Neutral allele theory considers polymor-

94

phism to be a transient phase of molecular
evolution (20). Under the assumptions of
this theory, the level of the between-species
divergence at different loci is positively cor­
related with the level of within-species poly­
morphism. This idea is embodied in a sim­
ple, two-by-two statistical test (HKA test)
(2 I) that calculates the expected amounts of
polymorphism and divergence at two loci
and can indicate whether the observed
amounts are consistent with the neutral
theory. Because we determined that there is
an excess of replacement substitutions, we
directed our attention to the silent variation
in jgw and Adh. Figure 4 shows the observed
numbers of silent substitutions and polymor­
phisms in the Adh-derived portion of jgw.
Whereas their within-species polymorphism
is comparable, the between-species diver­
gence of jgw is twofold greater than that of
Adh (x2 = 7.12, P < 0.01).

The Adh gene in these two species has a
very biased codon usage with a 84 to 87%
G+C content at the third codon positions.
Furthermore, the codon preference as shown
by the within-species polymorphism data in
D. yakuba (19) also shows a very high G+C
content (85. 7%). Nevertheless, the codon
usage in the Adh-derived portion of jgw
seems to be evolving to a smaller amount of
G+C. The G+C contents at the third
codon position summed over all silent poly­
morphism sites in the two species are 70.0%
for D. teissieri and 62.4% for D. yakuba,
respectively. The G+C content at the third
codon position in the fixed replacement sites
between the two species is 60.0%. This
evolution of a moderate G+C content of
third sites is consistent with the higher rate
of silent divergence for jgw and the general
observation of a negative correlation be­
tween the codon bias (and G+C content)
and the rate of silent divergence (22).

The results presented reveal the molecu­
lar, genetic, and evolutionary mechanisms

SCIENCE • VOL 260 • 2 APRIL 1993

that participated in the origin of the chimer­
ic processed functional gene jgw. The chi­
meric nature of the expressed product of this
novel gene is unique among processed genes,
including those few that are functional, such
as the second gene for rat insulin, the human
POK gene, and other examples that have
been interpreted as outcomes of retroposi­
tions (2, 23). It has been proposed that
retrotransposition can be a major cause of
intron loss during evolution (24), as it was
with jgw. The evolution of jgw also demon­
strates that retrotransposition can be a
source of new, intron-containing genes in
eukaryotic evolution. It is also unique
among instances of exon-shuffiing (25) in
that it clearly involved the retrotransposi­
tion of an mRNA. Experiments of cross­
hybridization with a genomic Southern
(DNA) blot and Northern analysis with the
use of a probe generated from the captured
portion of jgw indicated that the source of
the captured portion of jgw was a duplication
of an unrelated gene (26). The analysis of
the early molecular evolution of jgw indicat­
ed that jgw was functional from the begin­
ning and experienced strong adaptive evolu­
tion for what must have been a novel func­
tion. Prevailing theories of the origin of new
genes assume initial relaxation of selection
(I). Our results provide a contrary interpre­
tation in which natural selection is present
throughout the origin of a new gene.

REFERENCES AND NOTES

1 . S. Ohno, Evolution by Gene Duplication (Spring­
er-Verlag, Berlin, 1970); T. Ohta, Genome 31, 304
(1989); J. B. S. Haldane, The Causes of Evolution
(Longmans, New York, 1932).

2. W. Gilbert, Nature 271, 501 (1978); G. F. Hollis, P.
A. Hieter, 0. W. McBride, D. Swan, P. Leder, ibid.
296, 321 (1982); J. R. McCarrey and K. Thomas,
ibid. 326, 501 (1987).

3. In an ancient legend from China (San Hai Jing),
Jingwei, a daughter of the Emperor Yande (first
Chinese emperor 3000 B.C.), tragically drowned
while swimming in the East China Sea. Jingwei
was then reincarnated as a beautiful bird that
drops stones and wood into the sea in an attempt
to fill it, thus preventing others from drowning. We
used the name “jingwei” because this gene
avoided the usual fate of processed gene (death)
and was “reincarnated” into a new structure with
novel function.

4. C. H. Langley et al., Proc. Natl. Acad. Sci. U.S.A.
79, 5631 (1982).

5. P. Jeffs and M. Ashburner, Proc. R. Soc. London,
Ser. B 244, 151 (1991).

6. The sequences were completely determined from
both strands by the dideoxyribonucleotide chain­
terminating method [F. Sanger et al., Proc. Natl.
Acad. Sci. U.S.A. 74, 5463 (1977)] on single-str

PLoS BIOLOGY

Population Genomics: Whole-Genome Analysis
of Polymorphism and Divergence
in Drosophila simulans

1,2* 1,2* 1,2 3 1,2,4,5
David J. Begun , Alisha K. Holloway , Kristian Stevens , LaDeana W. Hillier , Yu-Ping Poh ,

6,7 6 8,9 1,2,10 11 12,13
Matthew W. Hahn , Phillip M. Nista , Corbin D. Jones , Andrew D. Kern , Colin N. Dewey , Lior Pachter ,

13 1,2*
Eugene Myers , Charles H. Langley

1 Department of Evolution and Ecology, University of California Davis, Davis, California, United States of America, 2 Center for Population Biology, University of California

Davis, Davis, California, United States of America, 3 Genome Sequencing Center, Washington University School of Medicine, St. Louis, Missouri, United States of America,

4 Institute of Molecular and Cellular Biology, National Tsing Hua University, Hsinchu, Taiwan Authority, 5 Research Center for Biodiversity, Academica Sinica, Taipei, Taiwan

Authority, 6 Department of Biology, Indiana University, Bloomington, Indiana, United States of America, 7 School of Informatics, Indiana University, Bloomington, Indiana,

United States of America, 8 Department of Biology, University of North Carolina, Chapel Hill, North Carolina, United States of America, 9 Carolina Center for Genome

Sciences, University of North Carolina, Chapel Hill, North Carolina, United States of America, 10 Center for Biomolecular Science and Engineering, University of California

Santa Cruz, Santa Cruz, California, United States of America, 11 Department of Biostatistics and Medical Informatics, University of Wisconsin, Madison, Wisconsin, United

States of America, 12 Department of Mathematics, University of California, Berkeley, California, United States of America, 13 Department of Computer Science, University of

California, Berkeley, California, United States of America

The population genetic perspective is that the processes shaping genomic variation can be revealed only through
simultaneous investigation of sequence polymorphism and divergence within and between closely related species.
Here we present a population genetic analysis of Drosophila simulans based on whole-genome shotgun sequencing of
multiple inbred lines and comparison of the resulting data to genome assemblies of the closely related species, D.
melanogaster and D. yakuba. We discovered previously unknown, large-scale fluctuations of polymorphism and
divergence along chromosome arms, and significantly less polymorphism and faster divergence on the X chromosome.
We generated a comprehensive list of functional elements in the D. simulans genome influenced by adaptive evolution.
Finally, we characterized genomic patterns of base composition for coding and noncoding sequence. These results
suggest several new hypotheses regarding the genetic and biological mechanisms controlling polymorphism and
divergence across the Drosophila genome, and provide a rich resource for the investigation of adaptive evolution and
functional variation in D. simulans.

Citation: Begun DJ, Holloway AK, Stevens K, Hillier LW, Poh YP, et al. (2007) Population genomics: whole-genome analysis of polymorphism and divergence in Drosophila
simulans. PLoS Biol 5(11): e310. doi:10.1371/journal.pbio.0050310

Introduction

Given the long history of Drosophila as a central model
system in evolutionary genetics beginning with the origins of
empirical population genetics in the 1930s, it is unsurprising
that Drosophila data have inspired the development of
methods to test population genetic theories using DNA
variation within and between closely related species [1–4].
These methods rest on the supposition of the neutral theory
of molecular evolution that polymorphism and divergence
are manifestations of mutation and genetic drift of neutral
variants at different time scales [5]. Under neutrality, poly-
morphism is a ‘‘snapshot’’ of variation, some of which
ultimately contributes to species divergence as a result of
fixation by genetic drift. Natural selection, however, may
cause functionally important variants to rapidly increase or
decrease in frequency, resulting in patterns of polymorphism
and divergence that deviate from neutral expectations [1,2,6].
A powerful aspect of inferring evolutionary mechanism in
this population genetic context is that selection on sequence
variants with miniscule fitness effects, which would be
difficult or impossible to study in nature or in the laboratory
but are evolutionarily important, may cause detectable
deviations from neutral predictions. Another notable aspect
of these population genetic approaches is that they facilitate

inferences about recent selection—which may be manifest as
reduced polymorphism or elevated linkage disequilibrium—
or about selection that has occurred in the distant past—
which may be manifest as unexpectedly high levels of
divergence. The application of these conceptual advances to
the study of variation in closely related species has resulted in
several fundamental advances in our understanding of the
relative contributions of mutation, genetic drift, recombina-
tion, and natural selection to sequence variation. However, it
is also clear that our genomic understanding of population
genetics has been hobbled by fragmentary and nonrandom
population genetic sampling of genomes. Thus, the full value

Academic Editor: Mohamed A. F. Noor, Duke University, United States of America

Received March 19, 2007; Accepted September 26, 2007; Published November 6,
2007

Copyright: � 2007 Begun et al. This is an open-access article distributed under the
terms of the Creative Commons Attribution License, which permits unrestricted
use, distribution, and reproduction in any medium, provided the original author
and source are credited.

Abbreviations: CDS, coding sequence; GO, gene ontology; indel, insertion/
deletion; MK test, McDonald and Kreitman test; UTR, untranslated region

* To whom correspondence should be addressed. E-mail: djbegun@ucdavis.edu
(DJB); akholloway@ucdavis.edu (AKH); chlangley@ucdavis.edu (CHL)

PLoS Biology | www.plosbiology.org 2534 November 2007 | Volume 5 | Issue 11 | e310

Population Genomics of D. simulans

Author Summary

Population genomics, the study of genome-wide patterns of
sequence variation within and between closely related species,
can provide a comprehensive view of the relative importance of
mutation, recombination, natural selection, and genetic drift in
evolution. It can also provide fundamental insights into the
biological attributes of organisms that are specifically shaped by
adaptive evolution. One approach for generating population
genomic datasets is to align DNA sequences from whole-genome
shotgun projects to a standard reference sequence. We used this
approach to carry out whole-genome analysis of polymorphism and
divergence in Drosophila simulans, a close relative of the model
system, D. melanogaster. We find that polymorphism and diver-
gence fluctuate on a large scale across the genome and that these
fluctuations are probably explained by natural selection rather than
by variation in mutation rates. Our analysis suggests that adaptive
protein evolution is common and is often related to biological
processes that may be associated with gene expression, chromo-
some biology, and reproduction. The approaches presented here
will have broad applicability to future analysis of population
genomic variation in other systems, including humans.

of genome annotation has not yet been applied to the study
of population genetic mechanisms.

Combining whole-genome studies of genetic variation
within and between closely related species (i.e., population
genomics) with high-quality genome annotation offers several
major advantages. For example, we have known for more than
a decade that regions of the genome experiencing reduced
crossing over in Drosophila tend to show reduced levels of
polymorphism yet normal levels of divergence between
species [7–10]. This pattern can only result from natural
selection reducing levels of polymorphism at linked neutral
sites, because it violates the neutral theory prediction of a
strong positive correlation between polymorphism and
divergence [5]. However, we have no general genomic
description of the physical scale of variation in polymor-
phism and divergence in Drosophila and how such variation
might be related to variation in mutation rates, recombina-
tion rates, gene density, natural selection, or other factors.
Similarly, although several Drosophila genes have been targets
of molecular population genetic analysis, in many cases, these
genes were not randomly chosen but were targeted because of
their putative association with phenotypes thought to have a
history of adaptive evolution [11,12]. Such biased data make it
difficult to estimate the proportion of proteins diverging
under adaptive evolution. In a similar vein, the unique power
of molecular population genetic analysis, when used in
concert with genome annotation, could fundamentally alter
our notions about phenotypic divergence due to natural
selection. This is because our current understanding of
phenotypic divergence and its causes is based on a small
and necessarily highly biased description of phenotypic
variation. Alternatively, a comprehensive genomic investiga-
tion of adaptive divergence could use genome annotations to
reveal large numbers of new biological processes previously
unsuspected of having diverged under selection. Here we
present a population genomic analysis of D. simulans. D.
simulans and D. melanogaster are closely related and split from
the outgroup species, D. yakuba, several million years ago [13–
15]. The vast majority of D. simulans and D. yakuba euchro-

matic DNA is readily aligned to D. melanogaster, which permits
direct use of D. melanogaster annotation for investigation of
polymorphism and divergence and allows reliable inference
of D. simulans–D. melanogaster ancestral states over much of the
genome. Our analysis uses a draft version of a D. yakuba
genome assembly (aligned to the D. melanogaster reference
sequence) and a set of light-coverage, whole-genome shotgun
data from multiple inbred lines of D. simulans, which were
syntenically aligned to the D. melanogaster reference sequence.

Results/Discussion

Genomes and Assemblies
Seven lines of D. simulans and one line of D. yakuba were

sequenced at the Washington University Genome Sequencing
Center (the white paper can be found at http://www.genome.
gov/11008080). The D. simulans lines were selected to capture
variation in populations from putatively ancestral geographic
regions [16], recent cosmopolitan populations, and strains
encompassing the three highly diverged mitochondrial
haplotypes previously described for the species [17]. These
strains have been deposited at the Tucson Drosophila Stock
Center (http://stockcenter.arl.arizona.edu). A total of 2,424,141
D. simulans traces and 2,245,197 D. yakuba traces from this
project have been deposited in the National Center for
Biotechnology Information (NCBI) trace archive. D. simulans
syntenic assemblies were created by aligning trimmed,
uniquely mapped sequence traces from each D. simulans strain
to the euchromatic D. melanogaster reference sequence (v4).
Two strains from the same population, sim4 and sim6, were
unintentionally mixed prior to library construction; reads
from these strains were combined to generate a single, deeper,
syntenic assembly (see Materials and Methods), which is
referred to as SIM4/6. The other strains investigated are
referred to as C167.4, MD106TS, MD199S, NC48S, and w501 .
Thus, six (rather than seven) D. simulans syntenic assemblies
are the objects of analysis. Details on the fly strains and
procedures used to create these assemblies, including the use
of sequence quality scores, can be found in Materials and
Methods. The coverages (in Mbp) for C167.4, MD106TS,
MD199S, NC48S, SIM4/6, and w501 , are 56.9, 56.3, 63.4, 42.6,
89.8, and 84.8, respectively. A D. yakuba strain Tai18E2 whole-
genome shotgun assembly (v2.0; http://genome.wustl.edu/)
generated by the Parallel Contig Assembly Program (PCAP)
[18] was aligned to the D. melanogaster reference sequence
(Materials and Methods). The main use of the D. yakuba
assembly was to infer states of the D. simulans–D. melanogaster
ancestor. For many analyses, we used divergence estimates for
the D. simulans lineage or the D. melanogaster lineage (from the
inferred D. simulans–D. melanogaster ancestor) rather than the
pairwise (i.e., unpolarized) divergence between these species.
These lineage-specific estimates are often referred to as ‘‘D.
simulans divergence,’’ ‘‘D. melanogaster divergence,’’ or ‘‘polar-
ized divergence.’’
A total of 393,951,345 D. simulans base pairs and

102,574,197 D. yakuba base pairs were syntenically aligned to
the D. melanogaster reference sequence. Several tens of
kilobases of repeat-rich sequences near the telomeres and
centromeres of each chromosome arm were excluded from
our analyses (Materials and Methods). D. simulans genes were
conservatively filtered for analysis based on conserved
physical organization and reading frame with respect to the

PLoS Biology | www.plosbiology.org 2535 November 2007 | Volume 5 | Issue 11 | e310

�

Population Genomics of D. simulans

Table 1. Autosome and X Chromosome Weighted Averages of Nucleotide Heterozygosity (p) and Lineage Divergence

Sequence Type Sites Chromosome p Divmel Divsim Divyak

Euchromatic Nonsynonymous X 0.0018 0.0067 0.0070 0.0253
A 0.0026 0.0061 0.0057 0.0223

Synonymous X 0.0199 0.0767 0.0519 0.2314
A 0.0352 0.0695 0.0524 0.2187

Intron X 0.0166 0.0248 0.0330 0.1175
A 0.0212 0.0240 0.0281 0.1028

59 UTR X 0.0079 0.0233 0.0258 0.1018
A 0.0108 0.0216 0.0203 0.0842

39 UTR X 0.0088 0.0199 0.0261 0.0957
A 0.0113 0.0186 0.0192 0.0775

Intergenic X 0.0153 0.0231 0.0299 0.1102
A 0.0204 0.0225 0.0265 0.0957

Heterochromatic Nonsynonymous X 0.0014 0.0088 0.0089 0.0269
A 0.0017 0.0083 0.0075 0.0354

Synonymous X 0.0132 0.0664 0.0493 0.2385
A 0.0136 0.0589 0.0523 0.2338

Divmel, D. melanogaster lineage divergence; Divsim, D. simulans lineage divergence; Divyak,
D. simulans/D. melanogaster common ancestor), see Materials and Methods.
doi:10.1371/journal.pbio.0050310.t001

D. melanogaster reference sequence gene models (Materials and
Methods). We took this conservative approach so as to retain
only the highest quality D. simulans data for most inferences.
The number of D. simulans genes remaining after filtering was
11,466. Ninety-eight percent of coding sequence (CDS)
nucleotides from this gene set are covered by at least one D.
simulans allele. The average number of lines sequenced per
aligned D. simulans base was 3.90. For several analyses in which
heterozygosity and divergence per site were estimated, we
further filtered the data so as to retain only genes or
functional elements (e.g., untranslated regions [UTRs]) for
which the total number of bases sequenced across all lines
exceeded an arbitrary threshold (see Materials and Methods).
The numbers of genes for which we estimated coding region
expected heterozygosity, unpolarized divergence, and polar-
ized divergence were 11,403, 11,439, and 10,150, respectively.
Coverage on the X chromosome was slightly lower than
autosomal coverage, which is consistent with less X chromo-
some DNA than autosomal DNA in mixed-sex DNA preps.
Variable coverage required analysis of individual coverage
classes (n ¼ 1–6) for a given region or feature, followed by
estimation and inference weighted by coverage (Materials and
Methods). The D. simulans syntenic alignments are available at
http://www.dpgp.org/. An alternative D. simulans ‘‘mosaic’’
assembly, which is available at http://www.genome.wustl.edu/,
was created independently of the D. melanogaster reference
sequence.

General Patterns of Polymorphism and Divergence
Nucleotide variation. We observed 2,965,987 polymorphic

nucleotides, of which 43,878 altered the amino acid sequence;
77% of sampled D. simulans genes were segregating at least
one amino acid polymorphism. The average, expected
nucleotide heterozygosity (hereafter, ‘‘heterozygosity’’ or
‘‘pnt ’’) for the X chromosome and autosomes was 0.0135 and
0.0180, respectively. X chromosome pnt was not significantly
different from that of the autosomes (after multiplying X
chromosome pnt by 4/3, to correct for X/autosome effective
population size differences when there are equal numbers of

D. yakuba lineage divergence (corresponds to divergence between D. yakuba and the

males and females; see [19]). However, X chromosome
divergence was greater than autosomal divergence in all
three lineages (50-kb windows; Table 1, Table S1, Figure 1,
Dataset S8). We will discuss this pattern in greater detail
below.
Not surprisingly, many patterns of molecular evolution

identified from previously published datasets were confirmed
in this genomic analysis. For example, synonymous sites and
nonsynonymous sites were the fastest and slowest evolving
sites types, respectively [20–24]. Nonsynonymous divergence
(dN) and synonymous divergence (dS) were positively, though
weakly, correlated (r2 ¼0.052, p , 0.0001) [25–27], and dN was
weakly, negatively correlated with CDS length (Spearman’s q
¼� 0.03, p ¼0.0005) [28,29]. More generally, longer functional
elements showed smaller D. simulans divergence than did
shorter elements (intron Spearman’s q ¼� 0.33; intergenic
Spearman’s q ¼� 0.39; 39 UTRs Spearman’s q ¼� 0.11: all show
p , 0.0001) [21,30].
Insertion/deletion (indel) variation. We investigated only

small indels (�10 bp), because they were inferred with high
confidence (Materials and Methods). Variants were classified
with respect to the D. melanogaster reference sequence;
divergence estimates were unpolarized. An analysis of trans-
posable element variation can be found in Text S1. Estimates
of small-indel heterozygosity for the X chromosome and
autosomes (Table S1) were lower than estimates of nucleotide
heterozygosity [31]. Interestingly, variation in nucleotide and
indel heterozygosity across chromosome arms was highly
correlated ([32], Figures 1 and 2; Spearman’s q between 0.45
and 0.69, p , 10 4 for each arm). Deletion heterozygosity and
divergence were consistently greater than insertion hetero-
zygosity and divergence (Figures S1 and S2, Datasets S11–S15)
for both the X chromosome and the autosomes, which
supports and extends previous claims, based on analysis of
repetitive sequences [33], of a general mutational bias for
deletions in Drosophila.
D. simulans autosomal pnt and divergence are of similar

magnitude. Mean polarized autosomal divergence (50-kb

PLoS Biology | www.plosbiology.org 2536 November 2007 | Volume 5 | Issue 11 | e310

Population Genomics of D. simulans

Figure 1. Patterns of Polymorphism and Divergence of Nucleotides along Chromosome Arms
Nucleotide p (blue) and div on the D. simulans lineage (red) in 150-kbp windows are plotted every 10 kbp. v[–log(p)] (olive) as a measure of deviation (þ
or –) in the proportion of polymorphic sites in 30-kbp windows is plotted every 10 kbp (see Materials and Methods). C and T correspond to locations of
centromeres and telomeres, respectively. Chromosome arm 3R coordinates correspond to D. simulans locations after accounting for fixed inversion on
the D. melanogaster lineage.
doi:10.1371/journal.pbio.0050310.g001

windows; 0.024) was only slightly greater than mean autoso-
mal pnt (0.018), even with regions of severely reduced pnt near
telomeres and centromeres included. Indeed, estimates of pnt
for several genomic regions are roughly equal to the genomic
average polarized divergence (Figure 1), suggesting the
existence of large numbers of shared polymorphisms in D.
simulans and D. melanogaster; such variants should be over-
represented in regions of higher nucleotide heterozygosity in
D. simulans. These patterns suggest that the average time to
the most recent common ancestor of D. simulans alleles is
nearly as great as the average time of the most recent
common ancestor of D. simulans and D. melanogaster. The
similarity in scale of polymorphism and divergence in D.
simulans also suggests that many of the neutral mutations that
have fixed in D. simulans were polymorphic in the common
ancestor of the two species. As we discuss below, this has
implications for interpreting chromosomal patterns of poly-
morphism and divergence in this species.

As expected under the neutral model, and given the
observation that much of the D. simulans lineage divergence
is attributable to polymorphism, D. simulans pnt and diver-
gence (50-kb windows) were highly, significantly correlated
(autosome Spearman’s q ¼ 0.56, p , 0.0001: X chromosome
Spearman’s q ¼ 0.48, p , 0.0001) [5]. Moreover, the genetic
and population genetic processes shaping patterns of
divergence along chromosome arms appear to operate in a
similar manner in D. simulans and D. melanogaster, as polarized
divergence (50-kb windows) for the two lineages was highly
correlated (Spearman’s q ¼ 0.74; p , 0.0001). Nevertheless,

some regions of the genome showed highly significant
increases in divergence in either the D. simulans or the D.
melanogaster lineage (see below).
Variation near centromeres and telomeres. Figure 1 and

Figure S1 support previous reports documenting severely
reduced levels of polymorphism in the most proximal and
distal euchromatic regions of Drosophila chromosome arms
[7,10,34–36]. The fact that divergence in such regions
(Materials and Methods) is only slightly lower (50-kb median
¼0.0238) than that of the rest of the euchromatic genome (50-
kb median ¼0.0248) (Mann-Whitney U, p , 0.0001), supports
the hypothesis that reduced pnt in these regions is due to
selection at linked sites rather than reduced neutral mutation
rates [1,3,6]. Genes that are located in repetitive regions of
chromosomes near telomeres and centromeres (Materials and
Methods), which we refer to as ‘‘heterochromatic,’’ showed
moderately reduced nonsynonymous and synonymous heter-
ozygosity compared with other genes (Table 1, Dataset S6)
[37] and showed a substantially higher ratio of nonsynon-
ymous-to-synonymous polymorphism and divergence relative
to other genes (Table S2) [38].
Interestingly, the magnitude and physical extent of reduced

pnt near telomeres and centromeres appears to vary among
arms. Moreover, the physical scale over which divergence
varied along the basal region of 3R appears to be much
smaller than the scale for other arms, which is seen in Figure
1 as a more compressed, thick red line representing
divergence. These heterogeneous patterns of sequence
variation near centromeres and telomeres across chromo-

PLoS Biology | www.plosbiology.org 2537 November 2007 | Volume 5 | Issue 11 | e310

Population Genomics of D. simulans

Figure 2. Patterns of Polymorphism for Nucleotides, Small Insertions, and Small Deletions along Chromosome Arms
p for nucleotides (blue), p for small (� 10 bp) insertions (orange), and p for small (� 10 bp) deletions (orchid) among the D. simulans lines in 150-kbp
windows are plotted every 10 kbp (see Materials and Methods). C and T correspond to locations of centromeres and telomeres, respectively.
Chromosome arm 3R coordinates correspond to D. simulans locations after accounting for fixed inversion on the D. melanogaster lineage.
doi:10.1371/journal.pbio.0050310.g002

some arms may reflect real differences. For example, genetic
data from D. melanogaster suggest that the centromere-
associated effects of reduced crossing-over are greater for
the autosomes than for the X chromosome and also suggest
that the X chromosome telomere is associated with a stronger
reduction in crossing-over compared with the autosomal
telomeres [39]. Alternatively, some of the heterogeneity
between chromosome arms in the centromere proximal
regions may reflect variation in the amount of repeat-rich
sequence excluded from the analysis (Materials and Methods).

X versus Autosome Divergence
Faster-X divergence. The X chromosome differs from the

autosomes in its genetics as well as in its population genetics
[40,41]. These differences have motivated several attempts to
compare patterns of polymorphism and divergence on these
two classes of chromosomes and to use such comparisons to
test theoretical population genetic models [19,41]. For
example, several population genetic models (e.g., recessivity
of beneficial mutations) predict faster evolution of X-linked
versus autosomal genes [42]. Nevertheless, there is currently
no statistical support for greater divergence of X-linked
versus autosomal genes in Drosophila [19,43,44].

The genomic data presented here clearly show that the X is
evolving faster than the autosomes. For example, median
(standard error [SE]) X versus autosome divergence for 50-kb
windows was 0.0274 (0.0003) versus 0.0242 (0.0001) for D.
simulans, 0.0233 (0.0002) versus 0.0223 (0.0007) for D. mela-
nogaster, and 0.1012 (0.0007) versus 0.0883 (0.0003) for D.
yakuba. The X evolves significantly faster than the autosomes in
D. simulans, D. melanogaster, and D. yakuba (Tables 1 and S1; 50-

kb windows, Mann-Whitney U; z ¼4.99, 12.92, and 14.68 for D.
melanogaster, D. simulans, and D. yakuba respectively, all p ,
0.0001), although the faster-X effect appeared to be consid-
erably smaller in D. melanogaster than in D. simulans or D. yakuba.
Moreover, of the 18 lineage divergence estimates (six site types
and three lineages), only one, D. simulans synonymous sites,
failed to show faster-X evolution (Table 1). However, not all
classes of site/lineages showed statistically significant faster-X
evolution (Table S3). Thus, the faster-X effect is likely to be
general for Drosophila but vary in magnitude across lineages
and site types. Mean X chromosome divergence in previous
analyses of smaller datasets [19,43,44] was higher (though not
significantly so) than autosome divergence, in agreement with
these genomic results. Finally, indel divergence also showed a
faster-X effect (Mann-Whitney U, p , 0.0001 for both
insertions and deletions).
Interestingly, the lengths of coding regions, introns, inter-

genic regions, and 59 and 39 UTRs were significantly longer
(Mann-Whitney U, all five have p , 0.0001) for the X
chromosome than for the autosomes in D. melanogaster [45].
Longer introns, intergenic sequences, and genes tend to
evolve more slowly than shorter functional elements (above
and [45]), suggesting that the faster-X inference is conserva-
tive. Perhaps the X chromosome requires additional sequen-
ces for proper regulation through dosage compensation (e.g.,
[46–48]) or proper large-scale organization in the nucleus
[49]. Alternatively, if directional selection were more com-
mon on the X chromosome, then Hill-Robertson effects [50]
could favor insertions, because selection is expected to be
more effective when there is more recombination between
selected sites. However, the fact that X-linked deletion

PLoS Biology | www.plosbiology.org 2538 November 2007 | Volume 5 | Issue 11 | e310

� �

� �

�

Population Genomics of D. simulans

divergence is much greater than insertion divergence, at least
for small indels (see below), does not support this idea.
Further analysis of larger indels could clarify this matter.
Finally, under the premise that ancestral polymorphism
makes a considerable contribution to D. simulans divergence,
lower X chromosome polymorphism (relative to ancestral
autosome polymorphism) would also make the faster-X
inference conservative.

As noted above, faster-X evolution has several possible
explanations, including recessivity of beneficial mutations,
underdominance, more frequent directional selection on
males than on females, higher mutation rates in females than
in males, or higher mutation rates on the X chromosome
versus the autosomes [19,40–42]. The fact that faster-X
evolution is observed across most site types is consistent with
the hypothesis that X chromosome mutation rates are greater
than autosomal mutation rates. The X chromosome is distinct
from the autosomes in that it is dosage compensated in males
through hypertranscription of X-linked genes [51–53]. Dosage
compensation of the Drosophila male germline [52] could
result in higher X-linked mutation rates if chromatin
conformation associated with hypertranscription increases
mutation rates. Indeed, cytological and biochemical studies of
the male Drosophila polytene chromosomes suggest that the X
has a fundamentally different chromatin organization than
the autosomes [54]. Alternatively, DNA repair in the hetero-
gametic male could have different properties than repair in
females. In addition to the possible contribution of elevated
X-linked mutation rates to faster-X evolution, some aspects of
the data support a role for selection in elevating X
chromosome substitution rates. For example, the three site
classes that showed the greatest X/autosome divergence ratio
in D. simulans (nonsynonymous, 59 UTR and 39 UTR) also
showed the strongest evidence for adaptive divergence in
contrasts of polymorphic and fixed variants in D. simulans (see
below). Furthermore, the observation of a significantly higher
frequency of derived polymorphic variants on the X relative
to the autosomes [55] (Table S4) is consistent with more
adaptive evolution on the X chromosome [56,57]. However,
there is no obvious enrichment of genes showing a history of
recurrent adaptive protein evolution on the X chromosome
(see below).

In addition to the overall faster rate of X chromosome
evolution, relative rate tests (Materials and Methods) revealed
that the deviations of observed numbers of substitutions from
neutral expectations are significantly greater for the X
chromosome than for autosomes in both D. simulans and D.
melanogaster (Mann-Whitney U, p ¼1.3 3 10 13 and 1.4 3 10 4
for D. simulans and D. melanogaster, respectively). The
magnitude of the deviations of D. simulans substitutions from
expected numbers (Materials and Methods) varied along
chromosome arms (Table S5 and Figure S3), with the X
chromosome showing a particularly strong physical clustering
of unusual regions. Though these patterns could be

Crystal Structure of an Ancient Protein: Evolution
by Conformational Epistasis
Eric A. Ortlund, et al.
Science 317, 1544 (2007);
DOI: 10.1126/science.1142819

The following resources related to this article are available online at
www.sciencemag.org (this information is current as of December 11, 2007 ):

Updated information and services, including high-resolution figures, can be found in the online
version of this article at:
http://www.sciencemag.org/cgi/content/full/317/5844/1544

Supporting Online Material can be found at:
http://www.sciencemag.org/cgi/content/full/1142819/DC1

A list of selected additional articles on the Science Web sites related to this article can be
found at:
http://www.sciencemag.org/cgi/content/full/317/5844/1544#related-content

This article cites 26 articles, 9 of which can be accessed for free:
http://www.sciencemag.org/cgi/content/full/317/5844/1544#otherarticles

This article has been cited by 1 article(s) on the ISI Web of Science.

This article has been cited by 1 articles hosted by HighWire Press; see:
http://www.sciencemag.org/cgi/content/full/317/5844/1544#otherarticles

This article appears in the following subject collections:
Evolution
http://www.sciencemag.org/cgi/collection/evolution

Information about obtaining reprints of this article or about obtaining permission to reproduce
this article in whole or in part can be found at:
http://www.sciencemag.org/about/permissions.dtl

D
o
w

n
lo

a
d
e
d
f
ro

m
w

w
w

.s
ci

e
n
ce

m
a
g
.o

rg
o

n
D

e
ce

m
b
e
r

1
1
,

2
0
0
7

Science (print ISSN 0036-8075; online ISSN 1095-9203) is published weekly, except the last week in December, by the
American Association for the Advancement of Science, 1200 New York Avenue NW, Washington, DC 20005. Copyright
2007 by the American Association for the Advancement of Science; all rights reserved. The title Science is a
registered trademark of AAAS.

REPORTS

proaches have been and are being considered. For
example, in Singapore, where 84% of the popu-
lation lives in public housing (35), regulations that
explicitly recognize the role of spatial segregation in
sectarianism specify the percentage of ethnic groups
to occupy housing blocks (36). This legally
compels ethnic mixing at a scale finer than that
which our study finds likely to lead to violence.
Given the natural tendency toward social separa-
tion, maintaining such mixing requires a level of
authoritarianism that might not be entertained in
other locations. Still, despite social tensions (37),
the current absence of violence provides some
support to our analysis. The alternative approach—
aiding in the separation process by establishing
clear boundaries between cultural groups to
prevent violence—has also gained recent atten-
tion (38, 39). Although further studies are
needed, there exist assessments (39) of the impact
of historical partitions in Ireland, Cyprus, the
Indian subcontinent, and the Middle East that
may be consistent with the understanding of type
separation and a critical scale of mixing or
separation presented here.

The insight provided by this study may help
inform policy debates by guiding our understanding
of the consequences of policy alternatives. The
purpose of this paper does not include promoting
specific policy options. Although our work re-
inforces suggestions to consider separation, we are
not diminishing the relevance of concerns about the
desirability of separation or its process. Even where
separation may be indicated as a way of preventing
violence, caution is warranted to ensure that the
goal of preventing violence does not become a
justification for violence. Moreover, even a peaceful
process of separation is likely to be objectionable.
There may be ways to positively motivate
separation using incentives, as well as to mitigate
negative aspects of separation that often include
displacement of populations and mobility barriers.

Our results for the range of filter diameters that
provide good statistical agreement between
reported and predicted violence in the former
Yugoslavia and India suggest that regions of width
less than 10 km or greater than 100 km may
provide sufficient mixing or isolation to reduce the
chance of violence. These bounds may be affected
by a variety of secondary factors including social
and economic conditions; the simulation resolu-
tion may limit the accuracy of the lower limit; and
boundaries such as rivers, other physical barriers,
or political divisions will surely play a role. Still,
this may provide initial guidance for strategic
planning. Identifying the nature of boundaries to
be established and the means for ensuring their
stability, however, must reflect local issues.

Our approach does not consider the relative
merits of cultures, individual acts, or immediate
causes of violence, but rather the conditions that may
promote violence. It is worth considering whether, in
places where cultural differentiation is taking place,
conflict might be prevented or minimized by political
acts that create appropriate boundaries suited to the
current geocultural regions rather than the existing

historically based state boundaries. Such bounda-
ries need not inhibit trade and commerce and need
not mark the boundaries of states, but should allow
each cultural group to adopt independent behav-
iors in separate domains. Peaceful coexistence
need not require complete integration.

References and Notes
1. M. White, Deaths by Mass Unpleasantness: Estimated

Total for the Entire 20th Century, http://users.erols.com/
mwhite28/warstat8.htm (September 2005).

2. D. L. Horowitz, Ethnic Groups in Conflict (Univ. of
California Press, Berkeley and Los Angeles, ed. 2, 2000).

3. B. Harff, T. R. Gurr, Ethnic Conflict in World Politics
(Westview, Boulder, ed. 2, 2004).

4. S. Huntington, The Clash of Civilizations and the Remaking of
World Order (Simon & Schuster, New York, 1996).

5. D. Chirot, M. E. P. Seligman, Eds., Ethnopolitical Warfare:
Causes, Consequences, and Possible Solutions (American
Psychological Association, Washington, DC, 2001).

6. M. Reynal-Querol, J. Conflict Resolut. 46, 29 (2002).
7. T. R. Gulden, Politics Life Sciences 21, 26 (2002).
8. H. Buhaug, S. Gates, J. Peace Res. 39, 417 (2002).
9. A. Varshney, Ethnic Conflict and Civic Life: Hindus and

Muslims in India (Yale Univ. Press, New Haven, CT, 2003).
10. M. D. Toft, The Geography of Ethnic Violence: Identity,

Interests, and the Indivisibility of Territory (Princeton
Univ. Press, Princeton, NJ, 2003).

11. J. Fox, Religion, Civilization, and Civil War: 1945 through the
New Millennium (Lexington Books, Lanham, MD, 2004).

12. M. Mann, The Dark Side of Democracy: Explaining Ethnic
Cleansing (Cambridge Univ. Press, New York, 2004).

13. I. Lustick, Am. Polit. Sci. Rev. 98, 209 (2004).
14. Materials and methods are available as supporting

material on Science Online.
15. T. C. Schelling, J. Math. Sociol. 1, 143 (1971).
16. J. Mimkes, J. Therm. Anal. 43, 521 (1995).
17. H. P. Young, Individual Strategy and Social Structure

(Princeton Univ. Press, Princeton, NJ, 1998).
18. R. Van Kempen, A. S. Ozuekren, Urban Stud. 35, 1631 (1998).
19. Y. Bar-Yam, in Dynamics of Complex Systems (Perseus

Press, Cambridge, MA, 1997), chap. 7.
20. A. J. Bray, Adv. Phys. 43, 357 (1994).
21. I. M. Lifshitz, V. V. Slyozov, J. Phys. Chem. Solids 19, 35

(1961).
22. D. A. Huse, Phys. Rev. B 34, 7845 (1986).

23. W. Easterly, R. Levine, Q. J. Econ. 112, 1203 (1997).
24. P. Collier, A. Hoeffler, Oxf. Econ. Pap. 50, 563 (1998).
25. R. H. Bates, Am. Econ. Rev. 90, 131 (2000).
26. J. D. Fearon, D. D. Laitin, Am. Pol. Sci. Rev. 97, 75 (2003).
27. D. N. Posner, Am. J. Pol. Sci. 48, 849 (2004).
28. I. Daubechies, Ten Lectures on Wavelets, (SIAM,

Philadelphia, 1992).
29. A. Arneodo, E. Bacry, P. V. Graves, J. F. Muzy, Phys. Rev.

Lett. 74, 3293 (1995).
30. P. Ch. Ivanov et al., Nature 383, 323 (1996).
31. Map of Yugoslavia, Courtesy of the University of Texas

Libraries., www.lib.utexas.edu/maps/europe/yugoslav.jpg.
32. R. Petrovic, Yugosl. Surv. 33, 3 (1992).
33. K. Chaudhuri, Frontline 18 (no. 2), www.hinduonnet.com/

fline/fl1802/18020330.htm.
34. Final Report, Carnegie Commission on Preventing Deadly

Conflict, www.wilsoncenter.org/subsites/ccpdc/pubs/
rept97/finfr.htm.

35. A. Brief Background, Housing and Development Board,
Singapore Government, www.hdb.gov.sg/fi10/fi10296p.nsf/
WPDis/About%20UsA%20Brief%20Background%20-%
20HDB’s%20Beginnings.

36. Ethnic Integration Policy, Housing and Development
Board, Singapore Government, www.hdb.gov.sg/fi10/
fi10201p.nsf/WPDis/Buying%20A%20Resale%
20FlatEthnic%20Group%20Eligibility.

37. D. Murphy, Christian Science Monitor, 5 February 2002,
www.csmonitor.com/2002/0205/p07s01-woap.html.

38. J. Tullberg, B. S. Tullberg, Politics Life Sciences 16, 237 (1997).
39. C. Kaufmann, Int. Secur. 23, 120 (1998).
40. We thank G. Wolfe, M. Woolsey, and L. Burlingame for

editing the manuscript; B. Wang for assistance with figures;
M. Nguyen and Z. Bar-Yam for assistance with identifying
data; and I. Epstein, S. Pimm, F. Schwartz, E. Downs, and
S. Frey for helpful comments. We acknowledge internal
support by the New England Complex Systems Institute and
the U.S. government for support of preliminary results.

Supporting Online Material
www.sciencemag.org/cgi/content/full/317/5844/1540/DC1
Methods
Figs. S1.1 to S4.3
SOM Text
Table S1
References
Bibliography

30 November 2006; accepted 13 August 2007
10.1126/science.1142734

Crystal Structure of an Ancient
Protein: Evolution by
Conformational Epistasis
Eric A. Ortlund,1* Jamie T. Bridgham,2* Matthew R. Redinbo,1 Joseph W. Thornton2†

The structural mechanisms by which proteins have evolved new functions are known only indirectly.
We report x-ray crystal structures of a resurrected ancestral protein—the ~450 million-year-old
precursor of vertebrate glucocorticoid (GR) and mineralocorticoid (MR) receptors. Using structural,
phylogenetic, and functional analysis, we identify the specific set of historical mutations that
recapitulate the evolution of GR’s hormone specificity from an MR-like ancestor. These
substitutions repositioned crucial residues to create new receptor-ligand and intraprotein contacts.
Strong epistatic interactions occur because one substitution changes the conformational position
of another site. “Permissive” mutations—substitutions of no immediate consequence, which
stabilize specific elements of the protein and allow it to tolerate subsequent function-switching
changes—played a major role in determining GR’s evolutionary trajectory.

D
o
w

n
lo

a
d
e
d
f
ro

m
w

w
w

.s
ci

e
n
ce

m
a
g
.o

rg
o

n
D

e
ce

m
b
e
r

1
1
,

2
0
0
7

A
central goal in molecular evolution is to
understand the mechanisms and dynam-
ics by which changes in gene sequence

generate shifts in function and therefore pheno-
type (1, 2). A complete understanding of this

process requires analysis of how changes in protein
structure mediate the effects of mutations on
function. Comparative analyses of extant proteins
have provided indirect insights into the diversifi-
cation of protein structure (3–6), and protein

1544 14 SEPTEMBER 2007 VOL 317 SCIENCE www.sciencemag.org

A B

F
o
ld

a
ct

iv
a
tio

n

30 4010HomoGR RajaGR HomoMR
8 3020
6

20
410

102
0 0
-10 -9 -8 -7 -6 -5 -11 -10 -9 -8 -7 -6 -5 -11 -10 -9 -8 -7 -6

0

Hormone (log M)

TetrapodGR TeleostGR ElasmobranchGR MRs(8)
(4) (6) (1) 20

AncGR2

~420 Ma

15

10

5

0
-11 -10 -9 -8 -7 -6

36aa
+1∆

20

15
AncGR1

10

25aa

~440 Ma
5

0
-11 -10 -9 -8

30
AncCR

20

-7 -6 C

C18

Aldosterone
Cortisol
DOC

10 C17
~470 Ma

0
-11 -10 -9 -8 -7 -6

C11

ormones.

REPORTS

engineering studies have elucidated structure-
function relations that shape the evolutionary
process (7–11). To directly identify the mecha-
nisms by which historical mutations generated
new functions, however, it is necessary to
compare proteins through evolutionary time.

Here we report the empirical structures of an
ancient protein, which we “resurrected” (12) by
phylogenetically determining its maximum likeli-
hood sequence from a large database of extant se-
quences, biochemically synthesizing a gene coding
for the inferred ancestral protein, expressing it in
cultured cells, and determining the protein’s
structure by x-ray crystallography. Specifically, we
investigated the mechanistic basis for the functional
evolution of the glucocorticoid receptor (GR), a
hormone-regulated transcription factor present in all
jawed vertebrates (13). GR and its sister gene, the
mineralocorticoid receptor (MR), descend from the
duplication of a single ancient gene, the ancestral
corticoid receptor (AncCR), deep in the vertebrate
lineage ~450 million years ago (Ma) (Fig. 1A) (13).
GR is activated by the adrenal steroid cortisol and
regulates stress response, glucose homeostasis, and
other functions (14). MR is activated by aldosterone
in tetrapods and by deoxycorticosterone (DOC) in
teleosts to control electrolyte homeostasis, kidney

1Department of Chemistry, University of North Carolina,
Chapel Hill, NC 27599, USA. 2Center for Ecology and
Evolutionary Biology, University of Oregon, Eugene, OR
97403, USA.

*These authors contributed equally to this work.
†To whom correspondence should be addressed. E-mail:
joet@uoregon.edu

Fig. 1. (A) Functional evolution

and colon function, and other processes (14). MR is
also sensitive to cortisol, though considerably less
so than to aldosterone and DOC (13, 15).
Previously, AncCR was resurrected and found to
have MR-like sensitivity to aldosterone, DOC, and
cortisol, indicating that GR’s cortisol specificity is
evolutionarily derived (13).

To identify the structural mechanisms by
which GR evolved this new function, we used
x-ray crystallography to determine the structures
of the resurrected AncCR ligand-binding domain
(LBD) in complex with aldosterone, DOC, and
cortisol (16) at 1.9, 2.0, and 2.4 Å resolution,
respectively (table S1). All structures adopt the
classic active conformation for nuclear receptors
(17), with unambiguous electron density for each
hormone (Fig. 1B and figs. S1 and S2). AncCR’s
structure is extremely similar to the human MR
[root mean square deviation (RMSD) = 0.9 Å for
all backbone atoms] and, to a lesser extent, to the
human GR (RMSD = 1.2 Å). The network of
hydrogen-bonds supporting activation in the
human MR (18) is present in AncCR, indicating
that MR’s structural mode of action has been
conserved for >400 million years (fig. S3).

Because aldosterone evolved only in the
tetrapods, tens of millions of years after AncCR,
that receptor’s sensitivity to aldosterone was
surprising (13). The AncCR-ligand structures
indicate that the receptor’s ancient response to
aldosterone was a structural by-product of its
sensitivity to DOC, the likely ancestral ligand,
which it binds almost identically (Fig. 1C). Key
contacts for binding DOC involve conserved

surfaces among the hormones, and no obligate
contacts are made with moieties at C11, C17, and
C18, the only variable positions among the three
hormones. These inferences are robust to uncer-
tainty in the sequence reconstruction: We modeled
each plausible alternate reconstruction [posterior
probability (PP) > 0.20] into the AncCR crystal
structures and found that none significantly af-
fected the backbone conformation or ligand inter-
actions. The receptor, therefore, had the structural
potential to be fortuitously activated by aldoster-
one when that hormone evolved tens of millions
of years later, providing the mechanism for evo-
lution of the MR-aldosterone partnership by mo-
lecular exploitation, as described (13).

To determine how GR’s preference for cortisol
evolved, we identified substitutions that occurred
during the same period as the shift in GR function.
We used maximum likelihood phylogenetics to de-
termine the sequences of ancestral receptors along
the GR lineage (16). The reconstructions had strong
support, with mean PP >0.93 and the vast majority
of sites with PP >0.90 (tables S2 and S3). We
synthesized a cDNA for each reconstructed LBD,
expressed it in cultured cells, and experimentally
characterized its hormone sensitivity in a reporter
gene transcription assay (16). GR from the com-
mon ancestor of all jawed vertebrates (AncGR1 in
Fig. 1A) retained AncCR’s sensitivity to aldoster-
one, DOC, and cortisol. At the next node, however,
GR from the common ancestor of bony vertebrates
(AncGR2) had a phenotype like that of modern
GRs, responding only to cortisol. This inference is
robust to reconstruction uncertainty: We introduced

D
o
w

n
lo

a
d
e
d
f
ro

m
w

w
w

.s
ci

e
n
ce

m
a
g
.o

rg
o

n
D

e
ce

m
b
e
r

1
1
,
2
0
0
7

of corticosteroid receptors. Dose-
response curves show transcrip-
tion of a luciferase reporter gene
by extant and resurrected ances-
tral receptors with varying doses
(in log M) of aldosterone (green),
DOC (orange), and cortisol (pur-
ple). Black box indicates evolution
of cortisol specificity. The number
of sequence changes on each
branch is shown (aa, replacement;
D, deletion). Scale bars, SEM of
three replicates. Node dates from
the fossil record (19, 20). For com-
plete phylogeny and sequences,
see fig. S10 and table S5. (B)
Crystal structure of the AncCR LBD
with bound aldosterone (green,
with red oxygens). Helices are la-
beled. (C) AncCR’s ligand-binding
pocket. Side chains (<4.2 Å from
bound ligand) are superimposed
from crystal structures of AncCR
with aldosterone (green), DOC
(orange), and cortisol (purple).
Oxygen and nitrogen atoms are
red and blue, respectively; dashed
lines indicate hydrogen bonds.
Arrows show C11, C17, and C18
positions, which differ among the h

www.sciencemag.org SCIENCE VOL 317 14 SEPTEMBER 2007 1545

20

15

10

AncGR1+
L111Q

AncGR1+
S106P, L111Q

0

5

10

15

20

5

0
-11 -10 -9 -8 -7 -6 -11 -10 -9 -8 -7 -6 -5

AncGR1+
S106P

0

5

10

15

20AncGR120

15

10

5

0
-11 -10 -9 -8 -7 -6 -11 -10 -9 -8 -7 -6 -5

REPORTS

plausible alternative states by mutagenesis, but
none changed function (fig. S4). GR’s specificity
therefore evolved during the interval between these
two speciation events, ~420 to 440 Ma (19, 20).

During this interval, there were 36 substitutions
and one single-codon deletion (figs. S5 and S6).
Four substitutions and the deletion are conserved in
one state in all GRs that descend from AncGR2 and
in another state in all receptors with the ancestral
function. Two of these—S106P and L111Q (21)—
were previously identified as increasing cortisol
specificity when introduced into AncCR (13). We
introduced these substitutions into AncGR1 and
found that they recapitulate a large portion of the
functional shift from AncGR1 to AncGR2, radi-
cally reducing aldosterone and DOC response
while maintaining moderate sensitivity to cortisol
(Fig. 2A); the concentrations required for half-
maximal activation (EC50) by aldosterone and
DOC increased by 169- and 57-fold, respectively,
whereas that for cortisol increased only twofold. A
strong epistatic interaction between substitutions
was apparent: L111Q alone had little effect on
sensitivity to any hormone, but S106P dramatically
reduced activation by all ligands. Only the
combination switched receptor preference from
aldosterone and DOC to cortisol. Introducing these
historical substitutions into the human MR yielded
a completely nonfunctional receptor, as did
reversing them in the human GR (fig. S7). These
results emphasize the importance of having the
ancestral sequence to reveal the functional impacts
of historical substitutions.

To determine the mechanism by which these
two substitutions shift function, we compared the
structures of AncGR1 and AncGR2, which were
generated by homology modeling and energy
minimization based on the AncCR and human
GR crystal structures, respectively (16). These
structures are robust to uncertainty in the recon-
struction: Modeling plausible alternate states did
not significantly alter backbone conformation,
interactions with ligand, or intraprotein interactions.
The major structural difference between AncGR1

Fig. 2. Mechanism for switching A
AncGR1’s ligand preference from al-

and AncGR2 involves helix 7 and the loop
preceding it, which contain S106P and L111Q
and form part of the ligand pocket (Fig. 2B and fig.
S8). In AncGR1 and AncCR, the loop’s position is
stabilized by a hydrogen bond between Ser106 and
the backbone carbonyl of Met103 . Replacing Ser106

with proline in the derived GRs breaks this bond
and introduces a sharp kink into the backbone,
which pulls the loop downward, repositioning and
partially unwinding helix 7. By destabilizing this
crucial region of the receptor, S106P impairs
activation by all ligands. The movement of helix
7, however, also dramatically repositions site 111,
bringing it close to the ligand. In this conforma-
tional background, L111Q generates a hydrogen
bond with cortisol’s C17-hydroxyl, stabilizing the
receptor-hormone complex. Aldosterone and DOC
lack this hydroxyl, so the new bond is cortisol-
specific. The net effect of these two substitu-
tions is to destabilize the receptor complex with
aldosterone or DOC and restore stability in a
cortisol-specific fashion, switching AncGR2’s pref-
erence to that hormone. We call this mode of
structural evolution conformational epistasis, be-
cause one substitution remodels the protein back-
bone and repositions a second site, changing the
functional effect of substitution at the latter.

Although S106P and L111Q (“group X” for
convenience) recapitulate the evolutionary switch
in preference from aldosterone to cortisol, the
receptor retains some sensitivity to MR’s ligands,
unlike AncGR2 and extant GRs. We hypothesized
that the other three strictly conserved changes that
occurred between AncGR1 and AncGR2 (L29M,
F98I, and deletion S212D) would complete the
functional switch. Surprisingly, introducing these
“group Y” changes into the AncGR1 and AncGR1 +
X backgrounds produced completely nonfunc-
tional receptors that cannot activate transcription,
even in the presence of high ligand concentrations
(Fig. 3A). Additional epistatic substitutions must
have modulated the effect of group Y, which pro-
vided a permissive background for their evolution
that was not yet present in AncGR1.

The AncCR crystal structure allowed us to
identify these permissive mutations by analyzing
the effects of group Y substitutions (Fig. 3B).
In all steroid receptors, transcriptional activity
depends on the stability of an activation-function
helix (AF-H), which is repositioned when the
ligand binds, generating the interface for tran-
scriptional coactivators. The stability of this
orientation is determined by a network of inter-
actions among three structural elements: the loop
preceding AF-H, the ligand, and helix 3 (17).
Group Y substitutions compromise activation be-
cause they disrupt this network. S212D eliminates
a hydrogen bond that directly stabilizes the AF-H
loop, and L29M on helix 3 creates a steric clash
and unfavorable interactions with the D-ring of
the hormone. F98I opens up space between helix
3, helix 7, and the ligand; the resulting instability
is transmitted indirectly to AF-H, impairing
activation by all ligands (Fig. 3B). If the protein
could tolerate group Y, however, the structures
predict that these mutations would enhance
cortisol specificity: L29M forms a hydrogen
bond with cortisol’s unique C17-hydroxyl, and
the additional space created by F98I relieves a
steric clash between the repositioned loop and
Met108 , stabilizing the key interaction between
Q111 and the C17-hydroxyl (Fig. 3B).

We hypothesized that historical substitutions
that added stability to the regions destabilized by
group Y might have permitted the evolving pro-
tein to tolerate group Y mutations and to complete
the GR phenotype. Structural analysis suggested
two candidates (group Z): N26T generates a new
hydrogen bond between helix 3 and the AF-H
loop, and Q105L allows helix 7 to pack more
tightly against helix 3, stabilizing the latter and,
indirectly, AF-H (Fig. 3B). As predicted, intro-
ducing group Z into the nonfunctional AncGR1 +
X + Y receptor restored transcriptional activity,
indicating that Z is permissive for Y (Fig. 3A).
Further, AncGR1 + X + Y + Z displays a fully
GR-like phenotype that is unresponsive to
aldosterone and DOC and maintains moderate

B

D
o
w

n
lo

a
d
e
d
f
ro

m
w

w
w

.s
ci

e
n
ce

m
a
g
.o

rg
o

n
D

e
ce

m
b
e
r

1
1
,
2
0
0
7

dosterone to cortisol. (A) Effect of
substitutions S106P and L111Q on the
resurrected AncGR1’s response to hor-
mones. Dashed lines indicate sensitivity

F
o
ld

a
ct

iv
a
tio

n

to aldosterone (green), cortisol (purple),
and DOC (orange) as the EC50 for
reporter gene activation. Green arrow
shows probable pathway through a
functional intermediate; red arrow,
intermediate with radically reduced
sensitivity to all hormones. (B) Struc-
tural change conferring new ligand
specificity. Backbones of helices 6 and
7 from AncGR1 (green) and AncGR2
(yellow) in complex with cortisol are
superimposed. Substitution S106P Hormone (log M)
induces a kink in the interhelical loop
of AncGR2, repositioning sites 106 and 111 (arrows). In this background, L111Q forms a new hydrogen bond with cortisol’s unique C17-hydroxyl (dotted red line).

1546 14 SEPTEMBER 2007 VOL 317 SCIENCE www.sciencemag.org

REPORTS

cortisol sensitivity. Both N26T and Q105L are
required for this effect (table S4). Strong epistasis
is again apparent: Adding group Z substitutions
in the absence of Y has little or no effect on ligand-
activated transcription, presumably because the
receptor has not yet been destabilized (Fig. 3A).
Evolutionary trajectories that pass through func-
tional intermediates are more likely than those
involving nonfunctional steps (22), so the only
historically likely pathways to AncGR2 are those
in which the permissive substitutions of group Z
and the large-effect mutations of group X occurred
before group Y was complete (Fig. 3C).

Fig. 3. Permissive substitutions in the
evolution of receptor specificity. (A)
Effects of various combinations of
historical substitutions on AncGR1’s
transcriptional activity and hormone-
sensitivity in a reporter gene assay.
Group Y (L29M, F98I, and S212D) abol-
ishes receptor activity unless groups X
(S106P, L111Q) and Z (N26T and
Q105L) are present; the XYZ combina-
tion yields complete cortisol-specificity.
The 95% confidence interval for each
EC50 is in parentheses. Dash, no acti-
vation. (B) Structural prediction of
permissive substitutions. Models of

AncGR1 (green) and AncGR2 (yellow)
are shown with cortisol. Group X and Y
substitutions (circles and rectangles)
yield new interactions with the C17-
hydroxyl of cortisol (purple) but de-
stabilize receptor regions required for
activation. Group Z (underlined) imparts
additional stability to the destabilized
regions. (C) Restricted evolutionary
paths through sequence space. The
corners of the cube represent states for
residue sets X, Y, and Z. Edges represent
pathways from the ancestral sequence
(AncGR1) to the cortisol-specific combi-

Our discovery of permissive substitutions in the
AncGR1-AncGR2 interval suggested that other
permissive mutations might have evolved even
earlier. We used the structures to pred

letters to nature

5. Charlesworth, B. The effect of background selection against deleterious mutations on weakly selected,
linked variants. Genet. Res. 63, 213±227 (1994).

6. Fay, J., Wycoff, G. J. & Wu, C.-I. Positive and negative selection on the human genome. Genetics 158,
1227±1234 (2001).

7. McDonald, J. H. & Kreitman, M. Adaptive evolution at the Adh locus in Drosophila. Nature 351, 652±
654 (1991).

8. Charlesworth, B., Morgan, M. T. & Charlesworth, D. The effect of deleterious mutations on neutral
molecular variation. Genetics 134, 1289±1303 (1993).

9. Maynard Smith, J. & Haigh, J. The hitch-hiking effect of a favourable gene. Genet. Res. 23, 23±35 (1974).
10. Begun, D. J. & Aquadro, C. F. levels of naturally occuring DNA polymorphism correlate with

recombination rates in D. melanogaster. Nature 356, 519±520 (1992).
11. Begun, D. The frequency distribution of nucleotide variation in Drosophila simulans. Mol. Biol. Evol.

18, 1343±1352 (2001).
12. Kliman, R. Recent selection on synonymous codon usage in Drosophila. J. Mol. Evol. 49, 343±351 (1999).
13. Adams, M. D. et al. The genome sequence of Drosophila melanogaster. Science 287, 2185±2195 (2000).
14. Powell, J. R. & DeSalle, R. Drosophila molecular phylogenies and their uses. Evol. Biol. 28, 87±138

(1995).
15. Haldane, J. B. S. The cost of natural selection. J. Genet. 55, 511±524 (1957).
16. Kimura, M. Evolutionary rate at the molecular level. Nature 217, 624±626 (1968).
17. Thompson, J. D., Higgins, D. G. & Gibson, T. J. ClustalWÐimproving the sensitivity of progressive

multiple alignment through sequence weighting, position-speci®c gap penalties and weight matrix
choice. Nucl. Acids Res. 22, 4673±4680 (1994).

18. Xia, X. Data Analysis in Molecular Biology and Evolution (Kluwer Academic, London, 2000).
19. Rozas, J. & Rozas, R. DnaSP version 3: an integrated program for molecular population genetics and

molecular evolution analysis. Bioinformatics 15, 174±175 (1999).
20. Yang, Z. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl.

Biosci. 13, 555±556 (1997).

Supplementary Information accompanies the paper on Nature’s website
(http://www.nature.com).

Acknowledgements
We thank B. Charlesworth, C.-I. Wu, S. Otto, M. Whitlock, T. Johnson, P. Awadalla,
J. Gillespie, G. McVean and P. Keightley for helpful discussions, and E. Moriyama for help
with data collection. N.G.C.S. was funded by the Biotechnology and Biological Sciences
Research Council (BBSRC) and A.E.-W. is funded by the Royal Society and the BBSRC.

Competing interests statement
The authors declare that they have no competing ®nancial interests.

Correspondence and requests for materials should be addressed to A.E.-W.
(e-mail: a.c.eyre-walker@sussex.ac.uk).

………………………………………………………..
Testing the neutral theory of
molecular evolution with
genomic data from Drosophila
Justin C. Fay*², Gerald J. Wyckoff*² & Chung-I Wu*³

* Committee on Genetics, University of Chicago, Chicago, Illinois 60637, USA
³ Department of Ecology and Evolution, University of Chicago, Chicago,
Illinois 60637, USA

…………………………………………………………………………………………………………………………….

Although positive selection has been detected in many genes, its
overall contribution to protein evolution is debatable1. If the bulk
of molecular evolution is neutral, then the ratio of amino-acid (A)
to synonymous (S) polymorphism should, on average, equal that
of divergence2. A comparison of the A/S ratio of polymorphism in
Drosophila melanogaster with that of divergence from Drosophila
simulans shows that the A/S ratio of divergence is twice as highÐa
difference that is often attributed to positive selection. But an
increase in selective constraint owing to an increase in effective
population size could also explain this observation, and, if so, all
genes should be affected similarly. Here we show that the differ-
ence between polymorphism and divergence is limited to only a

² Present addresses: Department of Genome Sciences, Lawrence Berkeley National Laboratory, Berkeley,
California 94720 (J.C.F.); Department of Human Genetics, University of Chicago, Chicago, Illinois 60637,
USA (G.J.W).

fraction of the genes, which are also evolving more rapidly, and this
implies that positive selection is responsible. A higher A/S ratio of
divergence than of polymorphism is also observed in other species,
which suggests a rate of adaptive evolution that is far higher than
permitted by the neutral theory of molecular evolution.
The neutral theory holds that the bulk of DNA divergence

between species is driven by mutation and drift, rather than by
positive darwinian selection3. But because the effect of positive
selection is often masked by negative selection4, detecting positive
selection is a challenging task. A rate of amino-acid substitution
greater than that of synonymous substitution can be explained only
by positive selection5, but such a criterion is very stringent as
negative selection lowers the rate of amino-acid substitution. A
high rate of amino-acid substitution is limited mostly to genes that
are involved in resistance to disease or in sexual reproduction, where
there is continual room for improvement6,7.
The McDonald±Kreitman test can detect positive selection even

in the presence of negative selection through a ratio of amino-acid
divergence to synonymous divergence greater than that of
polymorphism2. The A/S ratio of divergence is in¯ated above
polymorphism by advantageous amino-acid mutations, which
quickly sweep through a population but have a cumulative effect
on divergence. The McDonald±Kreitman test has been applied to
many genes individually, but only a few have yielded a signi®cant
excess of amino-acid divergence (Drosophila genes are reviewed in
refs 8, 9). This may in part be caused by a lack of power in detecting
positive selection in individual genes unless a large number of
adaptive substitutions have occurred.
For those genes that have yielded a signi®cant McDonald±

Kreitman test result, the A/S ratio of divergence is more than twice
as great as polymorphism10±12 . The effects of positive selection may
also be obscured by slightly deleterious amino-acid mutations
that in¯ate the A/S ratio of polymorphism but not divergence.
The effects of slightly deleterious mutations can be removed by
comparing common polymorphism with divergence, because dele-
terious amino-acid mutations are kept at low frequency in the
population4. This can only be done when the data from a large
number of genes are combined; individual genes rarely contain
more than a few common amino-acid polymorphisms.
An important but rarely appreciated assumption of the

McDonald±Kreitman test is that the selective constraint on a gene
remains constant over time. The selective constraint on a gene is
determined by the proportion of amino-acid mutations that are
deleterious3, 2Ns , -1, so both a change in the selection coef®cient
(s) and a change in effective population size (N) can result in a
change in selective constraint. Although it is well known that
selective constraint is not static across phylogenetic lineages13,14,
this assumption is rarely justi®ed in applications of the McDonald±
Kreitman test. Whereas the strength of selection on each gene might
¯uctuate over time depending on the genetic or environmental
background, a genome-wide change in constraint, such as that
caused by a change in effective population size, should produce a
consistent increase or decrease in the A/S ratio across all genes.
Alternatively, under positive selection each gene might be affected
to a different degree and some genes might not be affected at all.
To compare genomic patterns of amino-acid and synonymous

Table 1 Polymorphisms in D. melanogaster and divergence from D. simulans

Gene* Class Amino-acid Synonymous A/S
polymorphism, A polymorphism, S

………………………………………………………………………………………………………………………………………………………..
X-linked Rare (#12.5%) 4 67 0.06

Common (.12.5%) 6 46 0.13
Divergence 42 189 0.22

Autosomal Rare 79 126 0.63
Common 44 118 0.37
Divergence 421 521 0.81

………………………………………………………………………………………………………………………………………………………..
* There are 5 X-linked and 31 autosomal genes with a sample size of eight or greater (see text for the
data from all 45 genes).

1024 © 2002 Macmillan Magazines Ltd NATURE | VOL 415 | 28 FEBRUARY 2002 | www.nature.com

letters to nature

Table 2 African and non-African common polymorphism and divergence

Class Population Amino-acid Synonymous A/S
polymorphism, A polymorphism, S

………………………………………………………………………………………………………………………………………………………..
Polymorphism Non-African 48 124 0.39

African 40 159 0.25
Divergence 413 663 0.62
………………………………………………………………………………………………………………………………………………………..

site evolution, we tabulated polymorphism in D. melanogaster and
divergence from D. simulans from 45 gene surveys (Methods). If all
amino-acid and synonymous variation is neutral, then the A/S ratio
of polymorphism and divergence should be constant. The A/S ratio
of divergence (598/950 = 0.63) is signi®cantly greater than that of
common polymorphism (65/224 = 0.29; P , 10 -6). We compared
divergence with the common rather than the total polymorphism
because deleterious mutations at low frequency in¯ate the A/S ratio
of polymorphism. For the 36 genes with sample sizes of eight or
greater, there is a signi®cant excess of rare over common amino-acid
variation in autosomal genes (P = 0.022; Table 1), as is observed in
humans4. The absence of a difference in X-linked genes suggests that
the deleterious mutations are partially recessive and are more
readily eliminated from the X chromosome.
Both positive selection and an increase in selective constraint on

amino-acid changes can produce a higher A/S ratio of divergence
than of polymorphism. But only under certain restrictive conditions
is a genome-wide change in constraint possible. One such condition
is an increase in effective population size that is neither too distant
nor too recent in the evolutionary past. If this possibility can be
ruled out, positive selection may be the only viable explanation for
the high rate of amino-acid divergence.
If an increase in selective constraint resulted from a population

size increase associated with the spread of D. melanogaster outside
Africa15, it might be more appropriate to compare the A/S ratio of
the African population with that of divergence. Table 2, which
includes the 32 genes for which both African and non-African
populations were surveyed, shows that there is a signi®cantly larger
A/S ratio of divergence than of polymorphism in either population.
If a recent increase in effective population size increased constraint
on amino-acid polymorphism in both African and non-African
populations, then patterns of synonymous polymorphism might be
skewed towards rare variants. Neither African or non-African
populations show this pattern16. Finally, if there has been a decrease
in effective population size along the D. melanogaster lineage17,18, the
A/S ratio of polymorphism should be greater than that of divergence
between the two species.

12

10

8

6

4

2

0

–12 –8 –4 0 4 8 >10

ka < 0.02

ka > 0.02

Excess of amino-acid divergence

N
u
m

b
e
r

o
f

g
e
n
e
s

Figure 1 The distribution of the excess of amino-acid divergence contributed by each
gene. For reference, fast and slowly evolving genes are denoted by a rate of amino-acid
substitution (ka) greater than (®lled bars) or less than (open bars) 2%.

Table 3 Polymorphism and divergence in neutral and fast genes

Genes* Class Amino-acid Synonymous A/S
polymorphism, A polymorphism, S

………………………………………………………………………………………………………………………………………………………..
Neutral Rare 31 90 0.34

Common 16 69 0.23
Divergence 65 247 0.26

Fast Rare 48 36 1.33
Common 28 49 0.57
Divergence 356 274 1.30

………………………………………………………………………………………………………………………………………………………..
*X-linked genes are excluded.

If an increase in effective population size has produced a genome-
wide increase in selective constraint, the A/S ratio of all genes should
be affected. In Fig. 1, the distribution of each gene’s contribution to
the excess of amino-acid divergence suggests that there are two
classes of gene: neutral and rapidly evolving. The neutral class
comprises 34 genes that deviate by less than 10 amino-acid sub-
stitutions from that expected on the basis of the A/S ratio of all
common polymorphism. The remaining 11 genes all have a higher
A/S ratio of divergence than of polymorphism, and account for the
whole difference in the A/S ratio of polymorphism and divergence.
These genes are Acp26Aa, Acp29Ab, anon1A3, anon1E9, anon1G5, ci,
est-6, Ref2P, Rel, tra and Zw. As expected under positive selection,
which increases the rate of protein evolution, these 11 genes have a
high rate of amino-acid substitution (Fig. 1).
Can the pattern in Fig. 1 be explained by selection or demogra-

phy? Table 3 shows that, in the rapidly evolving genes, the A/S ratios
of divergence and of rare polymorphism are much higher than the
A/S ratio of the common polymorphism. This is expected if the
genes are under positive selection. Although a large increase in
population size in the recent past could account for the difference
between the A/S ratio of divergence and that of common poly-
morphism, this explanation is incompatible with the very small
difference found in the 26 neutral genes. Because both the neutral
and rapidly evolving genes have a higher A/S ratio of rare poly-
morphism than of common polymorphism, both should have been
affected by an increase in effective population size.
If positive selection is common, other species should also have an

A/S ratio of divergence greater than that of polymorphism. In
addition, any demographic scheme is not likely to be shared by
several species. In a study of eight genes in D. simulans, Drosophila
mauritiana and Drosophila sechellia, the A/S ratio of polymorphism
(A/S = 32/183) is 34% that of divergence (28/55)19. In a study of 42
genes with polymorphism in both D. melanogaster and D. simulans,
the A/S ratio of polymorphism is 65% that of divergence (N. G. C.
Smith and A. Eyre-Walker, personal communication). In another
study of 23 genes, the A/S ratio of polymorphism (45/305) is 30%
that of divergence along the D. simulans lineage (65/133)20. In
humans, the A/S ratio of common polymorphism (70/122) found
in 181 genes is 65% that of divergence (3,660/4,151) found in a
different set of 182 human and Old World monkey genes4.
Although these genomic patterns of variation are not explained

easily by the neutral theory, slightly deleterious mutations must
clearly be accounted for in attempting to measure positive selection.
In humans, 38% of amino-acid polymorphism was estimated to be
slightly deleterious4, and in D. melanogaster the estimate is 26%,
(0.63 – 0.37) ́ 126/123, from the combined neutral and rapidly
evolving genes (Table 3). These slightly deleterious mutations,
which are emphasized by the nearly neutral theory21, could
become effectively neutral and ®xed during a population bottleneck
of suf®cient severity, providing a burst of amino-acid substitutions
and an increase in the A/S ratio of divergence. We control for the
impact of these slightly deleterious mutations by comparing the
rapidly evolving class of gene to the neutral class (Fig. 1, Table 3).
Additional genomic data from other species will be needed to
estimate the general impact of these slightly deleterious mutations
on protein evolution. M

NATURE | VOL 415 | 28 FEBRUARY 2002 | www.nature.com © 2002 Macmillan Magazines Ltd 1025

letters to nature

Methods
Data
A literature search yielded 45 genes for which polymorphism had been surveyed in
D. melanogaster and for which an outgroup sequence was available. Of these, 36 had a
sample size of eight or greater, 32 had been surveyed in at least two African and two non-
African individuals and 10 were of X-linked genes. The 45 genes and their references are
listed in Supplementary Information.

Analysis
Polymorphism data was tabulated by hand or from GenBank accession numbers using
SITES21 or DNASP22. For each polymorphic site, the minor allele was classi®ed as rare
(# 12.5%) or common (. 12.5%). The cutoff of 12.5% was chosen to exclude deleterious
mutations from the common frequency class and to include those genes with samples of
eight or more in the analysis of rare compared to common polymorphism. Cutoffs of 10
and 15% produce similar results. We treated three alleles segregating at a single nucleotide
as two segregating sites and excluded complex variations. Divergence data was obtained by
comparing a randomly chosen sequence of D. melanogaster with that of D. simulans or, if
unavailable, either D. mauritiana or D. sechellia. The number of amino-acid and
synonymous substitutions between species was estimated using Kimura’s two-parameter
model to correct for multiple hits.

The contribution of each gene to the excess number of amino-acid substitutions was
calculated as the excess number of amino-acid substitutions minus the excess number of
amino-acid polymorphisms found in each gene. The excess for polymorphism and
divergence is A – S ́ (65/224), where A and S are the number of amino-acid and
synonymous substitutions, respectively, and 65/224 is the total number of amino-acid
polymorphisms divided by synonymous polymorphisms. (Ideally, the excess of amino-
acid divergence in each gene should be calculated using only polymorphism and
divergence in that gene but there is rarely suf®cient polymorphism in a single gene for
comparison with divergence.) We also calculated the contribution to the excess separately
for three groups of genes sorted by their rate of amino-acid divergence. The two methods
produced a similar distribution so the simpler method using a single group of genes was
used.

Received 27 June; accepted 4 December 2001.

1. Nei, M. Molecular Evolutionary Genetics (Columbia Univ. Press, New York, 1987).
2. McDonald, J. H. & Kreitman, M. Adaptive protein evolution at the Adh locus in Drosophila. Nature

351, 652±654 (1991).
3. Kimura, M. The Neutral Theory of Molecular Evolution (Cambridge Univ. Press, Cambridge, 1983).
4. Fay, J. C., Wyckoff, G. J. & Wu, C.-I. Positive and negative selection on the human genome. Genetics

158, 1227±1234 (2001).
5. Kimura, M. Preponderance of synonymous changes as evidence for the neutral theory of molecular

evolution. Nature 267, 275±276 (1977).
6. Yang, Z. & Bielawski, J. P. Statistical methods for detecting molecular adaptation. Trends Ecol. Evol. 15,

496±503 (2000).
7. Wyckoff, G. J., Wang, W. & Wu, C.-I. Rapid evolution of male reproductive genes in the descent of

man. Nature 403, 304±309 (2000).
8. Weinreich, D. M. & Rand, D. M. Contrasting patterns of nonneutral evolution in proteins encoded in

nuclear and mitochondrial genomes. Genetics 156, 385±399 (2000).
9. Moriyama, E. N. & Powell, J. R. Intraspeci®c nuclear DNA variation in Drosophila. Mol. Biol. Evol. 13,

261±277 (1996).
10. Eanes, W. F., Kirchner, M. & Yoon, J. Evidence for adaptive evolution of the G6pd gene in the

Drosophila melanogaster and Drosophila simulans lineages. Proc. Natl Acad. Sci. USA 90, 7475±7479
(1993).

11. Begun, D. J. & Whitley, P. Adaptive evolution of relish, a Drosophila NF-kB/IkB protein. Genetics 154,
1231±1238 (2000).

12. Tsaur, S. C., Ting, C. T. & Wu, C. I. Positive selection driving the evolution of a gene of male
reproduction, Acp26Aa, of Drosophila: II. Divergence versus polymorphism. Mol. Biol. Evol. 15, 1040±
1046 (1998).

13. Langley, C. H. & Fitch, W. M. An examination of the constancy of the rate of molecular evolution.
J. Mol. Evol. 3, 161±177 (1974).

14. Ohta, T. Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral
theory. J. Mol. Evol. 40, 56±63 (1995).

15. Lachaise, D. M., Cariou, M.-L., David, J. R., Lemeunier, F. & Tsacas, L. The origin and dispersal of the
Drosophila melanogaster subgroup: a speculative paleogeographic essay. Evol. Biol. 22, 159±225
(1988).

16. Andolfatto, P. Contrasting patterns of X-linked and autosomal nucleotide variation in Drosophila
melanogaster and Drosophila simulans. Mol. Biol. Evol. 18, 279±290 (2001).

17. Akashi, H. Codon bias evolution in Drosophila: Population genetics of mutation-selection drift. Gene
205, 269±278 (1997).

18. McVean, G. A., Vieira, J. Inferring parameters of mutation, selection and demography from patterns
of synonymous site evolution in Drosophila. Genetics 157, 245±257 (2001).

19. Kliman, R. M. et al. The population genetics of the origin and divergence of the Drosophila simulans
complex species. Genetics 156, 1913±1931 (2000).

20. Begun, D. J. The frequency distribution of nucleotide variation in Drosophila simulans. Mol. Biol. Evol.
18, 1343±1352 (2001).

21. Ohta, T. Slightly deleterious mutant substitutions during evolution. Nature 246, 96±98 (1973).
22. Hey, J. & Wakeley, J. A coalescent estimator of the population recombination rate. Genetics 145, 833±

846 (1997).
23. Rozas, J. & Rozas, R. DnaSP version 3: an integrated program for molecular population genetics and

molecular evolution analysis. Bioinformatics 15, 174±175 (1999).

Supplementary Information accompanies the paper on Nature’s website
(http://www.nature.com).

Acknowledgements
This work was supported by grants from the NIH and NSF to C.-I.W. and a Genetics
Training Grant and a Department of Education PhD fellowship to J.C.F.

Competing interests statement
The authors declare that they have no competing ®nancial interests.

Correspondence and requests for materials should be addressed to J.C.F.
(e-mail: jcfay@lbl.gov).

………………………………………………………..
Brain potential and functional MRI
evidence for how to handle two
languages with one brain
Antoni Rodriguez-Fornells*, Michael Rotte², Hans-Jochen Heinze²,
ToÈmme NoÈsselt² & Thomas F. MuÈ nte*

* Department of Neuropsychology, Otto von Guericke University,
UniversitaÈtsplatz 2, GebaÈude 24, 39106 Magdeburg, Germany
² Klinik fuÈr Neurologie 2, Otto von Guericke University, Leipzigerstrasse 44,
39120 Magdeburg, Germany

…………………………………………………………………………………………………………………………….

Bilingual individuals need effective mechanisms to prevent inter-
ference from one language while processing material in the other1.
Here we show, using event-related brain potentials and functional
magnetic resonance imaging (fMRI), that words from the non-
target language are rejected at an early stage before semantic
analysis in bilinguals. Bilingual Spanish/Catalan and monolingual
Spanish subjects were instructed to press a button when presented
with words in one language, while ignoring words in the other
language and pseudowords. The brain potentials of bilingual
subjects in response to words of the non-target language were
not sensitive to word frequency, indicating that the meaning of
non-target words was not accessed in bilinguals. The fMRI
activation patterns of bilinguals included a number of areas
previously implicated in phonological and pseudoword process-
ing2±5 , suggesting that bilinguals use an indirect phonological
access route to the lexicon of the target language to avoid
interference6.
High-pro®ciency bilingual subjects manage to understand and

speak one of their languages without apparent interference from the
other. This is a remarkable ability in the face of the fact that neuro-
imaging studies have revealed, at least for high-pro®ciency bilin-
guals, that neuro-anatomical representations of both languages are

Monolinguals

–2 m V

400 800 1,200 ms

Bilinguals

Spanish Catalan Pseudo

Figure 1 Lateralized readiness potentials (LRPs) from the main experiment indicating the
preparation of motor responses. The onset latency of the LRP to Spanish words, estimated
by the time at which the amplitude was signi®cantly different from zero for at least 4
consecutive time points (sequential t-tests)14, was 408 ms in the monolingual and 520 ms
in the bilingual group. No LRP activity is observed for Catalan words, indicating an
effective blocking of `word’ (go) responses in the bilingual group.

1026 © 2002 Macmillan Magazines Ltd NATURE | VOL 415 | 28 FEBRUARY 2002 | www.nature.com

Copyright © 1998 by the Genetics Society of America

Evidence for Genetic Hitchhiking Effect Associated With Insecticide
Resistance in Aedes aegypti

Guiyun Yan,* Dave D. Chadee† and David W. Severson*,1

* Department of Animal Health and Biomedical Sciences, University of Wisconsin, Madison, Wisconsin 53706 and
† Insect Vector Control Division, Ministry of Health, St. Joseph, Trinidad and Tobago, West Indies

Manuscript received April 3, 1997
Accepted for publication October 29, 1997

A B S T R A C T
Information on genetic variation within and between populations is critical for understanding the evo-

lutionary history of mosquito populations and disease epidemiology. Previous studies with Drosophila sug-
gest that genetic variation of selectively neutral loci in a large fraction of genome may be constrained by
fixation of advantageous mutations associated with hitchhiking effect. This study examined restriction frag-
ment length polymorphisms of four natural Aedes aegypti mosquito populations from Trinidad and Tobago,
at 16 loci. These populations have been subjected to organophosphate (OP) insecticide treatments for
more than two decades, while dichlor-diphenyltrichlor (DDT) was the insecticide of choice prior to this pe-
riod. We predicted that genes closely linked to the OP target loci would exhibit reduced genetic variation
as a result of the hitchhiking effect associated with intensive OP insecticide selection. We also predicted
that genetic variability of the genes conferring resistance to DDT and loci near the target site would be sim-
ilar to other unlinked loci. As predicted, reduced genetic variation was found for loci in the general chro-
mosomal region of a putative OP target site, and these loci generally exhibited larger FST values than other
random loci. In contrast, the gene conferring resistance to DDT and its linked loci show polymorphisms
and genetic differentiation similar to other random loci. The reduced genetic variability and apparent
gene deletion in some regions of chromosome 1 likely reflect the hitchhiking effect associated with OP in-
secticide selection.

MOSQUITOES are important vectors for several human pathogens because of their close associa-
tion with humans. Mosquito habitats often change rap-
idly as a result of vector control efforts; therefore,
successful adaptation to varying human habitats is es-
sential for mosquito reproduction. Adaptation ability
of an organism depends on its genetic variability. Infor-
mation on genetic variation within and between popu-
lations is critical for understanding the evolutionary
history of mosquito populations and disease epidemiol-
ogy (Tabachnick and Black 1996). Protein electro-
phoresis and DNA sequence analyses have revealed
remarkable variation in many genes in natural popula-
tions of Drosophila and other species, but the genetic
variability seems to differ substantially for genes in dif-
ferent genome regions (Aquadro 1992). Distribution
patterns of genetic variants in natural populations are
the joint effects of various evolutionary forces and de-
mographic factors, including random genetic drift, se-
lection, recombination, mutation, gene flow, mating
system and life history (e.g., colonization, range expan-
sions or contractions; Slatkin 1985). Population life

Corresponding author: Guiyun Yan, Department of Biological Sci-
ences, State University of New York, 109 Cooke Hall, Buffalo, NY
14260. E-mail: gyan@calshp.cals.wisc.edu

1Present address: Department of Biological Sciences, State Univer-
sity of Notre Dame, Notre Dame, IN 46556.

history and mating structure influence all loci equally,
but selection affects only the target loci (Kreitman and
Akashi 1995). Variation of selectively neutral loci may
also be constrained by the hitchhiking effect, particu-
larly in genome regions with low recombination rates
and under extensive selection (Maynard Smith and
Haigh 1974). Recent studies with D. melanogaster sug-
gest that the hitchhiking effect may have occurred over
a large fraction of the insect genome (Begun and
Aquadro 1992). In this study we analyzed restriction
site variation of 16 loci in natural populations of the
yellow fever mosquitoes, Aedes aegypti, and provide evi-
dence that the hitchhiking effect may have reduced ge-
netic variation in the genome regions around a puta-
tive insecticide resistance locus.

A. aegypti is an important vector of yellow fever and
dengue fever viruses in many tropical countries, includ-
ing Trinidad and Tobago, West Indies. Control efforts
for A. aegypti have focused primarily on habitat reduc-
tion and chemical treatment, which is based on the de-
struction of breeding sites and the use of insecticides,
including dichlor-diphenyltrichlor (DDT) in the 1950s
and several organophosphates (OP) since the 1960s.
The wide use of insecticides has been a powerful selec-
tion agent, and rapid development of resistance to
DDT and OPs is well documented (Gilkes et al. 1956;
Rawlins and Wan 1995). The genetic mechanisms of
insect resistance to various insecticides have been well

Genetics 148: 793–800 (February, 1998)

794 G. Yan et al.

characterized. For example, a point mutation in the
para sodium channel gene confers one form of resis-
tance to DDT (Williamson et al. 1996), and esterase
(EST ) gene amplification is associated with resistance
to OPs in Culex mosquitoes (Mouchès et al. 1990). A
genetic linkage map, based largely on random cDNA
sequences, has been constructed for A. aegypti (Sever-
son et al. 1993), and several insecticide resistance genes
have been mapped (Severson et al. 1997). In this study,
we used cDNA markers distributed across the mosquito
genome to examine DNA polymorphism and popula-
tion genetic differentiation, and to examine mosquito
genome structural changes associated with strong selec-
tion imposed by insecticides. Several studies have
investigated the genetic variation of various genera of
mosquito populations with isozyme, RAPD-DNA, mic-
rosatellite and mitochondrial DNA markers (Powell et
al. 1980; Tabachnick and Wallis 1985; Conn et al.
1993; Chevillon et al. 1995; Apostol et al. 1996). Re-
striction fragment length polymorphism (RFLP) mark-
ers are particularly suitable for population genetic stud-
ies, because they are presumably neutral, highly
polymorphic, segregate as codominant markers, and
can be used for studies of other mosquito species (Sev-
erson et al. 1994a). We chose Trinidad and Tobago
populations because the population history is known
and surveillance programs have been well established
there. Population historical information is important
for the interpretation of genetic data. Because the mos-
quito populations have been under selection of OP in-
secticides, genetic variation of loci closely linked to an
esterase locus conferring resistance would be reduced
if the hitchhiking effect has occurred. The hitchhiking
effect would be prominent in genome regions with
strong linkage disequilibrium and intense selection
(e.g., insecticides). In contrast, gene diversity at the
para locus and other neighboring loci is expected to be
similar to unlinked loci in the genome, because selec-
tion pressure has been removed at the para locus since
DDT was abandoned more than two decades ago.

M AT E R I A L S A N D M E T H O D S

Natural history of A. aegypti in Trinidad and Tobago: It is gen-
erally believed that domestic A. aegypti originated from an
African sylvan ancestor, and was introduced to the New World
from West Africa via transoceanic trade during the fifteenth
to seventeenth centuries (Tabachnick 1991). Caribbean pop-
ulations probably represent the initial introduction of the
mosquito species into the New World in the course of New
World colonization. The first outbreak of yellow fever in Trin-
idad was recorded in 1796, and in 1820 in Tobago.

In the 1950s, intensive vector control programs aimed to-
ward mosquito eradication were adopted, primarily by the
widespread usage of DDT. In the early 1960s, Trinidad was
considered free of A. aegypti, but was reinfested in 1962. Over-
all, mosquito populations in Tobago have been exposed to
fewer insecticides than Trinidad populations. A. aegypti is now
widely distributed in Trinidad and Tobago, despite continued

intensive vector control efforts through the use of OP insecti-
cides. Insecticide applications not only impose strong selec-
tion on the target loci, but also lead to recurrent reductions
of population sizes.

Collection of samples: In conjunction with the A. aegypti
surveillance program, in April, 1995, we collected three geo-
graphically-distinct samples from Trinidad and one sample
from Tobago (Figure 1). These four villages share similar cli-
mates, including temperature and the annual amount of rain-
fall. For each village, 100 ovitraps were distributed, and about
half of the village’s residential area was covered (approxi-
mately two traps every five houses). Each ovitrap consisted of
a black plastic container roughly half-filled with water into
which a rectangular masonite strip was placed in an upright
manner. Female A. aegypti mosquitoes will readily oviposit on
the masonite strip, near the water interface. After 2–3 days,
the masonite strips were removed, transported to the labora-
tory, where attached eggs were allowed to hatch, and reared
into adults. All adults were identified as A. aegypti by micro-
scopic examination, and were frozen for subsequent DNA
analysis. Previous studies with A. aegypti in Puerto Rico suggest
that the mean number of families represented per ovitrap was
4.7 (e.g., several female mosquitoes frequently oviposit in the
same container; Apostol et al. 1994). Therefore, it is unlikely
that siblings within a subpopulation would be sampled.

RFLP and probe selection: We genotyped a total of 870
mosquitoes for four populations (n 5 150 for Curepe, 262 for
Couva, 258 for San Fernando, and 200 for Tobago). DNA ex-
traction from individual mosquitoes, digestion with EcoRI,
Southern blotting and hybridization were as previously de-
scribed (Severson et al. 1993). Fifteen mapped RFLP markers
were selected to provide broad coverage of the A. aegypti ge-
nome with an average resolution of 10.6 cM (Figure 2). All
clones used were random cDNA clones with the exception of
the para and Mal I clones. Mal I is a gene specifically expressed
in the salivary glands, and its putative function is related to

Figure 1.—Map of Trinidad and Tobago. Three samples
were collected from Trinidad: Curepe (10838.62’N,
61824.23’W), Couva (10826.12’N, 61828.19’W), and San
Fernando (10818.11’N, 61828.21’W). One sample was col-
lected from Tobago (11811.23’N, 60844.21’W).

795 Hitchhiking Effect in Mosquitoes

Figure 2.—Relative map positions of the 16 Aedes aegypti
(2N 5 6) RFLP loci used in the study. Chromosome numbers
are in italics. Map distances are in Kosambi centimorgans.
Underlined loci were not used in the study. Esterase gene
amplification is involved with resistance to organophosphate
insecticides in mosquitoes (Mouchès et al. 1990). Three
esterase loci were mapped to the general chromosomal loca-
tions as shown in the figure (Munstermann 1990). LF250
represents duplicated loci. Markers in italics are genes with
known functions; other markers are random cDNA.

sugar metabolism ( James et al. 1989). The mosquito genome
consists of single or low-copy DNA sequences and repetitive
DNA with short-period interspersion (Black and Rai 1988).
In this study, we focused on allelic variations of single- or low-
copy cDNA sequences.

Data analysis: DNA polymorphism and Hardy-Weinberg equilib-
rium (HWE) tests: Molecular weights of fragments detected by
each clone were estimated by comparing them to lambda-
HindIII digest standards included on each gel, using the Eagle
Sight image capture and analysis software (Stratagene, La
Jolla, CA). DNA polymorphisms may be measured by the
proportion of polymorphic loci, number of alleles, and het-
erozygosity. Conformance with HWE was tested using the
probability test for each locus and each population, using
the GENEPOP computer program (Raymond and Rousset
1995). Because this test is robust to allele frequencies, rare al-
leles were not pooled. We further tested whether distortion
from HWE resulted from deficient or excessive heterozygos-
ity, using the FIS statistics (Weir 1990; Rousset and Raymond
1995). FIS is defined as [1 2 (observed heterozygosity/ex-
pected heterozygosity from HWE)]. Because FIS estimates at
individual loci may be unduly influenced by rare alleles, we
tested the significance of the average FIS over all loci using the
method of Robertson and Hill (1984). Variations in het-
erozygosity among the populations were analyzed following
the method of Weir (1990), using the analysis of variance
(ANOVA) with subpopulations, individuals, loci and inter-
actions of loci and individuals as factors. All factors were
treated as random effects except loci.

Population genetic structure, gene flow and genetic distance:
Population genetic structure was examined with Wright’s
F-statistics, based on the procedure of Weir and Cockerham
(1984) and using the FSTAT computer program (Goudet
1995). Standard deviations (SD) of F-statistics were obtained

for each locus by a jackknife procedure over the alleles, and
were used to test the significance of the F statistics. We first
tested whether the three populations from Trinidad were sig-
nificantly substructured, then included the Tobago popula-
tion data in the analysis.

Gene flow (Nm) was estimated from the standardized-
among-population genetic variance (FST) estimate of each lo-
cus using the relationship Nm 5 (1/FST 2 1)/4, where N is the
effective population size of a deme, and m is the rate of gene
flow (Wright 1943). This equation assumes the infinite-is-
land model of population structure and gene flow. Few popu-
lations probably conform to this assumption, but it provides a
useful approximation of the relative magnitude of gene flow.
Gene flow was also estimated using the private-alleles method
for the appropriate loci (Slatkin and Barton 1989). Private
alleles are the alleles unique to a given deme. Nei’s unbaised
genetic distance for all pairs of populations was calculated
based on population allele frequency for all loci (Nei 1987).

R E S U LT S

DNA polymorphisms and HWE tests: Fifteen cDNA
markers examined in this study were all polymorphic.
The RFLP patterns of one marker (LF250) indicate
that this marker represents a gene duplication (data
not shown), and therefore, the 15 markers represented
a total of 16 loci. A total of 91 unique alleles were iden-
tified, 68 alleles (74.7%) were common to all four pop-
ulations. The average number of alleles was about five
per locus (Table 1). Six loci (LF198, ARC1, LF250a,
para, LF168 and Mal I) exhibited private alleles, and
five private alleles were in the Tobago population. An
excess of rare alleles was found: 16 alleles (18%) had a
frequency less than 0.05. Under the infinite alleles
model (equation 8.24; Kimura 1983), we expected to
find only four alleles in this frequency class with our
sample size (n 5 870). The mean sizes of restriction
fragments detected by the cDNA clones weighted by
their frequencies ranged from 0.73 kb at locus LF250b to
12.90 kb at locus LF352, and exhibited an overall mean
of 5.13 kb (95% confidence interval: 3.61–6.64 kb).

In general, high heterozygosity was observed in all
four mosquito populations, except at the LF90 locus. The
LF90 locus showed significantly lower heterozygosity
than the other 15 loci examined (Table 1; ANOVA, t 5
8.37, d.f. 5 1, P , 0.0001). The most heterozygous loci
were LF178 on chromosome 1 and LF282 on chromo-
some 2. The high heterozygosity at the LF178 locus does
not seem to be a result of sex linkage (see Figure 2), be-
cause males and females showed similar heterozygosity
(data not shown). Population average heterozygosity over
all 16 loci varied little among populations (ranged from
0.582 for Couva to 0.627 for Tobago), and such variations
were not statistically significant (Table 1; ANOVA, F 5
1.08, d.f. 5 3, 49, P . 0.05). Heterozygosity is not corre-
lated with the mean size of restriction fragments
weighted by frequencies at a locus (r 5 0.22, d.f. 5 15,
P . 0.05), but seems to correlate with the number of ob-
served alleles (r 5 0.49, d.f. 515, P 5 0.052).

The genotype frequencies at several loci did not con-

796 G. Yan et al.

TABLE 1

RFLP polymorphisms of four Aedes aegypti populations from Trinidad and Tobago, measured by
observed heterozygosity and the number of alleles

Curepe Couva San Fernando Tobago

Chromosome Locus n Hobs
aFIS n Hobs FIS n Hobs FIS n Hobs FIS

1 LF90 5 0.122 20.036 3 0.159 20.076 5 0.256 0.160* 3 0.391 0.244***
LF230b 3 0.250 0.600*** 3 0.069 0.880*** 3 0.129 0.774*** 3 0.296 0.535***
LF198 6 0.655 0.096 6 0.504 0.100** 6 0.775 0.029 7 0.738 20.120**
LF178 6 0.761 0.012 6 0.849 20.075* 6 0.851 20.100*** 6 0.828 20.120**
TY7 5 0.503 0.160* 5 0.546 0.151* 5 0.543 20.001 5 0.656 0.060*

Average over
chromosome 1 5.5 0.510 0.059 5.0 0.515 0.025 5.5 0.606 0.022 5.3 0.654 0.016

2 ARC1 5 0.750 20.046 5 0.724 0.015 5 0.702 20.038 6 0.725 0.010
LF138 4 0.729 20.083 4 0.714 20.201** 4 0.492 20.075 4 0.487 0.145**
LF282 9 0.838 20.009 7 0.795 0.025 8 0.864 20.060* 7 0.793 20.102*
LF98 5 0.642 0.092 5 0.654 20.035 6 0.682 20.073 6 0.878 20.068*
LF250a 4 0.615 0.060 3 0.596 20.032 4 0.800 20.154* 3 0.562 20.183***
LF250b 3 0.644 20.009 3 0.687 20.038 3 0.662 20.238*** 3 0.557 20.160*
LF115 5 0.514 20.031 5 0.366 20.016 5 0.566 20.101 4 0.562 20.077

Average over
chromosome 2 5.0 0.676 20.003 4.6 0.648 20.040 5.0 0.681 20.105 4.7 0.652 20.062

3 LF352 6 0.548 0.268** 6 0.727 0.113** 6 0.441 0.292*** 6 0.550 0.263***
LF261 4 0.483 0.035 4 0.279 20.061 4 0.476 20.058 3 0.592 0.008
para 4 0.469 0.152 4 0.563 0.024 5 0.598 0.081 6 0.555 0.103
LF168 6 0.667 0.042** 7 0.516 0.145** 7 0.667 0.109 6 0.582 20.043
MalI 4 0.627 0.103** 4 0.637 0.009 3 0.640 0.034 4 0.577 0.137*

Average over
chromosome 3 4.8 0.559 0.119 5.0 0.544 0.046 5.0 0.564 0.092 5.0 0.571 0.093

Average over all loci 5.1 0.598 0.050 4.8 0.582 0.003 5.1 0.626 20.012 4.9 0.627 0.006

*P , 0.05, ** P , 0.01, ***P , 0.001. Hobs, observed heterozygosity; n, number of alleles.
a Significant FIS also indicates distortion from HWE. Positive FIS indicates heterozygosity deficit from HWE expectation; negative

FIS indicates excess of heterozygosity.
b A large proportion of individuals show an apparent gene deletion at the LF230 locus, therefore, this locus was not used for

chromosomal average heterozygosity calculation. The heterozygosity and FIS were based on the individuals without deletions. The
percentage of individuals showing the gene deletion at the LF230 locus was 41.4 for the Curepe population, 58.6 for Couva, 53.7
for San Fernando, and 46.0 for Tobago.

form to HWE. Loci on chromosome 2 generally exhib-
ited a heterozygote excess, but loci on chromosome 3
that showed HWE distortion exhibited a heterozygote
deficit (Table 1). The FIS values varied greatly among
the loci, suggesting no systematic inbreeding occurred
in these populations. The average FIS over all loci was
not significantly different from 0 for each population
(Table 1). Departure from HWE probably reflects
either the effect of insecticide selection on some loci
linked to the resistance loci, or simply sampling error.

Population genetic structure, gene flow and genetic
distance: Analysis of F statistics for the three Trinidad
populations found small, but statistically significant FST
estimates for all loci (Table 2), suggesting that these
populations are genetically differentiated. FST estimates
showed a six-fold difference among loci, with an aver-
age FST over all loci of 0.043. When the Tobago popula-
tion was included in the analysis, the basic pattern of

estimation among the loci was not altered (Table FST

2). As expected, slightly larger FST values were obtained
for most loci, and the average FST over all loci was
0.056. The para locus exhibited similar polymorphism
and genetic differentiation as other random loci.

Assuming that the populations are at an equilibrium
between migration and random drift, the average num-
ber of migrants exchanged per generation can be cal-
culated. Average gene flow (Nm) among the four popu-
lations, based on the FST method, was 4.2 migrants per
generation (95% confidence interval: 3.2–5.7). This es-
timate was similar to the estimate based on the average
frequency of six private alleles present in the popula-
tions (Nm 5 4.5). Table 3 shows genetic distances and
gene flow between each pair of populations calculated
from the pair-wise average FST. A large gene flow be-
tween the Tobago and Trinidad populations was de-
tected. There was no significant correlation between ge-
netic distance and geographic distance (r2 5 0.5, d.f. 5
5, P . 0.05).

797 Hitchhiking Effect in Mosquitoes

TABLE 2

FST statistics and Nm estimates of four Aedes aegypti populations of Trinidad and Tobago

FST Nm estimates of all populations
a

Chromosome Locus
Trinidad

populationsb
Trinidad and Tobago

populationsa
Based on

FST
Based on the

private-alleles method

1

2

3

LF90
LF198
LF178
TY7
ARC1
LF138
LF282
LF98
LF250a
LF250b
LF115
LF352
LF261
para
LF168
MalI

0.053 6 0.022
0.121 6 0.062
0.011 6 0.006
0.033 6 0.022
0.064 6 0.046
0.044 6 0.032
0.019 6 0.011
0.042 6 0.012
0.081 6 0.053
0.106 6 0.067
0.021 6 0.012
0.119 6 0.059
0.038 6 0.034
0.039 6 0.025
0.020 6 0.013
0.034 6 0.019

0.107 6 0.102
0.109 6 0.046
0.025 6 0.013
0.034 6 0.011
0.049 6 0.034
0.030 6 0.022
0.041 6 0.027
0.061 6 0.021
0.099 6 0.045
0.110 6 0.041
0.055 6 0.039
0.079 6 0.037
0.102 6 0.060
0.035 6 0.017
0.040 6 0.021
0.020 6 0.015

2.1
2.0
9.8
7.1
4.9
8.1
5.9
3.9
2.3
2.0
4.3
2.9
2.2
6.9
6.0

12.3

—c
3.5

—
—

177.6
—
—
—
0.8

—
—
—
—
70.9
1.7
9.2

Summary over 16 loci 0.043 6 0.007 0.056 6 0.008 4.2 4.5

Values are 6 SD.
All FST values were significantly larger than 0 at P , 0.001. The test was performed using a jackknifing proce-

dure over samples.
a n 5 4.
b n 5 3.
c The estimate was not available because no private alleles existed for the locus.

Hitchhiking effect on DNA polymorphisms: Gene-
tic hitchhiking occurs when a (neutral) mutation
changes frequency through genetic linkage to a muta-
tion that is selected, resulting in reduced genetic varia-
tion surrounding the target site of selection. Low DNA
polymorphism at the LF90 locus suggests that hitchhik-
ing has probably occurred in the genome region of the
EST-4 locus. We collected additional evidence to test for
this hypothesis by examining genetic polymorphism at
the LF230 locus, which also is in the general genomic
region of EST-4 (Figure 2). An apparent chromosomal
deletion event occurred around this locus in 42–59%
of the individuals (Figure 3), and low heterozygosity
(0.07–0.29) was observed among those individuals with-
out the apparent deletion (Table 1). However, substan-
tial reduction of heterozygosity for the para locus and
other loci in the vicinity of para was not observed (Ta-
ble 1).

D I S C U S S I O N

DNA polymorphisms of four A. aegypti mosquito
populations were examined using RFLP markers. The
Trinidad populations have been exposed to OPs every
3–4 months for about two decades. These populations
have therefore experienced intense selection by insecti-
cides, that probably resulted in periodic population
bottlenecks. A population bottleneck maintains a long-

term effect on population heterozygosity, even for spe-
cies with a large intrinsic rate of growth such as A. aegypti
(Nei et al. 1975). Thus, genetic polymorphisms are ex-
pected to decline rapidly during insecticide use for any
locus in the mosquito genome. Loci conferring OP re-
sistance are expected to maintain lower genetic vari-
ability than other random loci in the genome, and ge-
netic variation of other closely-linked neutral loci may
be reduced through genetic linkage.

If recurrent population bottlenecks have occurred
in the mosquito populations, low polymorphism for all
loci in the genome would be expected. In contrast to
the expectation, we observed high polymorphisms for
most loci. For the same loci, average heterozygosities of

TABLE 3

Nm matrix based on pairwise FST estimates and Nei’s
unbiased genetic distance matrix

Curepe Couva San Fernando Tobago

Curepe
Couva
San Fernando
Tobago

11.4
8.6
4.3

0.041

3.7
2.6

0.057
0.112

3.7

0.113
0.166
0.124

Numbers below diagonal line are FST estimates, above the
diagonal line are Nei’s unbiased genetic distance.

798 G. Yan et al.

Figure 3.—Southern blot analysis of natural Aedes aegypti
populations probed with cDNA clones LF230 (top) and LF90
(bottom). The mosquito genomic DNA was digested with
EcoRI. Each lane is for a single mosquito. (Top) Apparent
gene deletion around the LF230 locus in 55% of individuals
(11 out of 20). (Bottom) Probe LF90 was used as a control to
demonstrate that absence of hybridization of mosquito ge-
nomic DNA to LF230 was not due to incomplete DNA diges-
tion, or to poor probe conditions. DNA hybridization was ob-
served for the same individuals with all other markers tested.
See Table 1 for heterozygosity and percentage of gene dele-
tion at the LF230 locus.

the populations studied here are substantially higher
than a laboratory population, which has not been ex-
posed to insecticides for more than 20 years and has
not experienced population bottleneck (Yan et al.
1997). The highest heterozygosity was observed for loci
at two chromosomal regions (LF198-LF178 on chromo-
some 1, and LF282-LF98 on chromosome 2). Coinci-
dentally, these two chromosomal regions in A. aegypti
harbor genes determining vector competence for filar-
ial worms and malaria parasites (Severson et al. 1994b,
1995). The high levels of heterozygosity observed may
be explained by two mechanisms. The first is that the
effective size of population bottlenecks has never been
small, because heterogeneous habitats may provide
effective shelters for the field populations. The sec-
ond is that genetic polymorphisms are introduced and
maintained by large gene flow among populations. The
gene flow estimates seem to support the second hy-
pothesis.

The LF90 locus consistently exhibited lower het-
erozygosity than other loci in the genome for the four
populations used in this study. The heterozygosity of a
RFLP locus may be influenced by several factors, in-
cluding the size of the probes, the size of the regions
being probed by the probes, reduced mutation or re-
combination rates in these genome regions, natural or
artificial selection on a particular locus, and hitchhik-
ing (selective sweep) of a selectively neutral locus by se-
lectively favored substitutions at linked loci. We argue

that the polymorphism pattern of the LF90 locus likely
reflects the result of a hitchhiking effect. First, the puta-
tive function of the LF90 clone is coding for ribosomal
protein S14 (Severson and Zhang 1996), and thus the
RFLP fragments of LF90 themselves are presumably
neutral. Second, LF90 is located in the general chro-
mosomal region of EST-4, and gene amplification at an
esterase locus is a common mechanism of OP resis-
tance. Our populations have been under selective pres-
sure by OP insecticides for decades. Third, low het-
erozygosity of the LF90 locus is likely not related to the
size of the genome region being probed by the LF90
marker, because we found no significant correlation be-
tween heterozygosity and RFLP fragment sizes. Fourth,
the observed heterozygosity of the LF90 locus in labora-
tory colonies of A. aegypti that have not been exposed to
insecticide selection was similar to other random loci
across the genome (Yan et al. 1997).

Our argument for a hitchhiking effect is strength-
ened by the RFLP data of the LF230 locus, which is
linked to LF90 and also is in the general genomic re-
gion of EST-4. We found very low DNA polymorphism
and an apparent gene deletion for many individuals at
this locus, a phenomenon which has not been observed
in other laboratory colonies of A. aegypti (Yan et al.
1997). Gene deletions may be the result of unequal re-
combination in this chromosomal region, associated
with the esterase gene amplification. For example, OP
resistant Culex mosquitoes possess 250–500 copies of a
30-kb esterase B1 gene, compared to a single copy in
susceptible individuals (Mouchès et al. 1990). Given
the genome size of A. aegypti of about 320 Mb (Zaitlin
and Severson, unpublished results), if the magnitude
of esterase gene amplification in the Aedes mosquitoes
is similar to Culex mosquitoes, then the hitchhiking ef-
fect associated with OP insecticide selection could af-
fect meiotic pairing across a large genome region of
chromosome 1, and could lead to chromosomal abnor-
malities (i.e., deletions or duplications) within this re-
gion. In contrast, a hitchhiking effect associated with
DDT usage is not evident for the para locus, as indi-
cated by the fact that genetic heterozygosity at the para
locus and loci closely linked to para was similar to other
random loci in the genome. This result is consistent
with the hypothesis that in the years since DDT was aban-
doned, the populations have had time to re-equilibrate.

Ideally, the hitchhiking effect should be demon-
strated at the nucleotide diversity level (Kaplan et al.
1989). Unfortunately, nucleotide diversity cannot be
appropriately estimated for the present data, because
our RFLP data is based on one restriction enzyme. To
statistically rule out the possibility of low heterozygosity
caused by reduced mutation rates and/or increased
functional constraints in the LF90 gene region, one
needs to examine intraspecific variation and interspe-
cific divergence over several gene regions for closely-
related species. This method has been elegantly applied

799 Hitchhiking Effect in Mosquitoes

to Drosophila studies (Begun and Aquadro 1991,
1992). The rationale is that, if reduced mutation rate in
a gene region leads to low intraspecific variation, then
interspecific divergence is expected to be smaller than
in other gene regions. However, hitchhiking effect only
redu

Open access, freely available online PLoS BIOLOGY

Two Rounds of Whole Genome Duplication
in the Ancestral Vertebrate

1 1,2
Paramvir Dehal , Jeffrey L. Boore

1 Evolutionary Genomics Department, Department of Energy Joint Genome Institute and Lawrence Berkeley National Laboratory, Walnut Creek, California, United States of

America, 2 Department of Integrative Biology, University of California, Berkeley, California, United States of America

The hypothesis that the relatively large and complex vertebrate genome was created by two ancient, whole genome
duplications has been hotly debated, but remains unresolved. We reconstructed the evolutionary relationships of all
gene families from the complete gene sets of a tunicate, fish, mouse, and human, and then determined when each
gene duplicated relative to the evolutionary tree of the organisms. We confirmed the results of earlier studies that
there remains little signal of these events in numbers of duplicated genes, gene tree topology, or the number of genes
per multigene family. However, when we plotted the genomic map positions of only the subset of paralogous genes
that were duplicated prior to the fish–tetrapod split, their global physical organization provides unmistakable
evidence of two distinct genome duplication events early in vertebrate evolution indicated by clear patterns of four-
way paralogous regions covering a large part of the human genome. Our results highlight the potential for these large-
scale genomic events to have driven the evolutionary success of the vertebrate lineage.

Citation: Dehal P, Boore JL (2005) Two rounds of whole genome duplication in the ancestral vertebrate. PLoS Biol 3(10): e314.

Introduction

It has long been hypothesized that the increased complex-
ity and genome size of vertebrates has resulted from two
rounds (2R) of whole genome duplication (WGD) occurring in
early vertebrate evolution, thus providing the requisite raw
materials [1]. This seemed to be supported by the long-
standing speculation that humans have about 100,000 genes,
roughly four times the number expected for invertebrates’
genomes, but this is now known to be incorrect, with the
actual human gene count being closer to 30,000 [2,3].
Conflicting analyses have now made this very controversial,
with some studies supporting 2R (e.g., [4–8]), others seeing
only a single round of WGD (e.g., [9–11]), and still others
refuting WGD altogether by concluding that nothing greater
than limited segmental duplications have occurred (e.g.,
[12,13]).

The 2R hypothesis had been bolstered by observations that
a few gene families, e.g., Hox clusters [14], follow a ‘‘4:1 rule’’
in the numbers of vertebrate to invertebrate genes. However,
comparison of the complete genome sequences of human
[2,3] and Drosophila [15] revealed that less than 5% of
homologous gene families follow the 4:1 rule [12]. Further,
although two sequential duplications are expected to
generate the evolutionary topology (AB)(CD) for the descend-
ent genes, rather than (A)(BCD), in fact, the relationships of
vertebrate multigene families do not generally show this
pattern, as indicated by early studies using only a few genes
[16] and confirmed as complete genome sequences became
available [2,13]. (However, for a different view, see [17].)
Several studies have incorporated data from sparse sampling
of genes from taxa thought to have branched near to these
purported duplications, including lamprey [18], hagfish,
amphioxus [17,19–22], and Ciona [23]; although these results
are useful for timing duplications, the conclusions could
never be viewed as definitively resolving this issue because
these products could have alternatively been generated by
duplications of individual genes or short gene segments

rather than by WGDs. Even duplicating all of the genes in a
genome individually is quite different from a whole genome
duplicating simultaneously.
There are several reasons why this has been a difficult issue

to resolve. After duplication, only the minority of gene pairs
will adopt a new function (‘‘neofunctionalization’’) or
partition old functions (‘‘subfunctionalization’’) quickly
enough to escape disabling mutations that would lead to
their eradication [24]; therefore, rampant gene loss rapidly
erases this signal of genome duplication. Further, four-
member gene families, even those with the (AB)(CD) top-
ology, can be generated by two rounds of duplications of
individual genes or of segments much smaller than the entire
genome, generating a condition that cannot be differentiated
on this basis from 2R followed by many gene losses. This
alternative scenario seems especially plausible because recent
analyses have shown that gene duplications occur much more
frequently than had been thought, with the typical rate being
sufficient to duplicate an entire genome equivalent every 100
million years (MY) [25,26]. Until recently, no complete
genome sequence has been available from an outgroup that
is closely related to vertebrates, and all methods of
phylogenetic reconstructions are less accurate with more
distant relatives such as Drosophila and yeast [20]. Lastly, there
has not been to date a method to accurately and compre-
hensively cluster genes into homologous families because
methods that rely on sequence similarity alone are highly

Received September 7, 2004; Accepted July 8, 2005; Published September 6, 2005
DOI: 10.1371/journal.pbio.0030314

Copyright: � 2005 Dehal and Boore. This is an open-access article distributed
under the terms of the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided the
original work is properly cited.

Abbreviations: 2R, two rounds of whole genome duplication; MY, million years;
WGD, whole genome duplication

Academic Editor: Peter Holland, University of Oxford, United Kingdom

E-mail: PSDehal@LBL.gov (PD), JLBoore@LBL.gov (JLB)

PLoS Biology | www.plosbiology.org 1700 October 2005 | Volume 3 | Issue 10 | e314

Vertebrate Genome Duplication

subject to artifactual association of slowly evolving paralogs
and to erroneous exclusion of the more rapidly evolving
genes.

Fortunately, as has been shown convincingly for the yeast
genome and for Arabidopsis [27–30], evidence of an ancient
genome duplication can be seen in the large-scale pattern of
the physical locations of homologous genes, even when the
great majority of the duplicated genes have been lost. Studies
have shown that the human genome also has multiple regions
of colinear paralogous gene copies [4,21,22,31–37], but
considered the arrangements of too small a number of genes
and genomic regions to be comprehensive. This approach is
now available for a large-scale evaluation of the vertebrate 2R
hypothesis, because complete (i.e., at least draft quality)
genome sequences are available for the tunicate Ciona
intestinalis [38] (a basal chordate outgroup) and the verte-
brates Takifugu rubripes [39] (a pufferfish ‘‘fugu’’), Mus musculus
[40] (mouse), and human [2,3]. Figure 1 illustrates how the

Figure 1. Pattern Predicted for the Relative Locations of Paralogous
Genes from Two Genome Duplications
(A) Representation of a hypothetical genome that has 22 genes shown as
colored squares.
(B) A genome duplication generates a complete set of paralogs in
identical order.
(C) Many paralogous genes suffer disabling mutations, become
pseudogenes, and are then lost. One could imagine this condition being
evidence of a single round of genome duplication followed by significant
gene losses.
(D) A second genome duplication recreates another set of paralogs in
identical order, with multigene families that retained two copies now
present in four, and those that had lost a member now present in two
copies.
(E) Again, many paralogous genes suffer disabling mutations, become
pseudogenes, and are then lost. Of course, unrelated gene duplications
and transpositions can occur. Even though this leaves only a few four-
member gene families, the patterns of 2- and 3-fold gene families unite,
in various combinations, all four genomic segments, revealing that the
sequential duplications had been of very large regions, in this case all or
nearly all of this hypothetical genome.
DOI: 10.1371/journal.pbio.0030314.g001

signal of two rounds of genome duplication could be retained
by the large-scale pattern in location of duplicated genes, in
which many tracks of paralogous duplicates (which may not
contain identical subsets of genes) each occur at exactly four
positions in the genome, i.e., ‘‘tetra-paralogons.’’ No similar
signal would be generated by repeated duplications of genes
or even large gene segments; only WGDs would result in such
global organization of paralogous genes.

Results

Gene Clustering and Duplication Timing
A graph-based method was used with the complete gene

sets of the four chordates (98,517 total genes; see Table 1 for
details of each step in the analysis) to generate clusters such
that each contains all, and only, those genes that descended
from a single gene in their common ancestor (Figure 2). A
multiple sequence alignment and a maximum likelihood
evolutionary tree were constructed for each cluster, then a
Web browser interface was built so that each can be viewed
individually. (For more details and updates that include more
taxa, see the ‘‘PhIGs’’ [Phylogenetically Inferred Groups] Web
site at http://phigs.org/.) We could then easily determine when
each gene duplicated relative to lineage splitting by compar-
ing these gene trees with the known evolutionary relation-
ships of the animals. For example, a gene duplication that is
specific to only one animal’s lineage is seen as two genes from
the same genome clustered together. A gene that duplicated
once in the unique common ancestor of mouse and human
would generate a tree that groups gene copy 1 of human and
mouse and, separately, gene copy 2 of human and mouse. Put
more generally, gene duplications that are shared by more
that one species are seen as a replication of the phylogeny of
the descendant organisms for each gene copy. Of course, gene
losses and various combinations of these processes are seen as
well. Figure 3 shows all possible gene topologies along with
how each would be interpreted.

This reveals that 46.6% of the ancestral chordate genes
appear in duplicate in one or more of the vertebrate lineages,
with 34.5% having at least one duplication before the
divergence of fish from tetrapods and 23.5% having at least
one duplication afterward. (Some of these are counted twice,
having had duplications both before and after the fish–
tetrapod split.) This means that there are 3,753 gene
duplications placed at the base of Vertebrata, which is
remarkable because the ancestral genome would be reason-
ably estimated to have had fewer than 20,000 genes, which is
the case for the tunicate as well as other invertebrate
outgroups. However, as can be seen in Figure 4, gene
duplications are in large numbers on every branch of the
tree, making it unclear whether this, in itself, indicates a
significant acceleration in duplication rate. Additionally, of
the gene clusters with duplications basal to the fish–tetrapod
split, 20.6% have had one duplication event, 10.8% have had
two, and 5.1% have had more than two, counter to the
expectation from 2R, and casting further doubt on the
significance of this for the 2R hypothesis.

Gene Family Membership
An early observation in support of 2R was that several gene

families have expanded from a single member in inverte-
brates to having four members for some vertebrates. Previous

PLoS Biology | www.plosbiology.org 1701 October 2005 | Volume 3 | Issue 10 | e314

Vertebrate Genome Duplication

Table 1. Overview of the Process for Analyzing the Complete Gene Sets with Number of Genes Included at Each Step

Step Process Step Gene Counts Clusters
Number

Ciona Fugu Mouse Human

1 Retrieve sequences 15,852 37,241 22,444 22,980 —
2 Run BLASTP 12,448 27,090 20,918 20,718 —
3 Make seeds — — — — —
4 Generate clusters 7,438 11,339 10,069 10,290 6,641
5 With duplication in vertebrates 3,623 8,394 7,131 7,235 3,096
6 With at least one fugu and one 3,402 7,885 6,907 7,015 2,951

tetrapod gene, ,100 copies in
each taxon

7 Create multiple sequence alignment — — — — —
8 Trim gaps, alignability 2,565 5,618 4,924 4,987 2,340
9 Perform phylogenetic analysis — — — — —

10 Determine strictly bifurcating nodes 1,776 3,770 3,118 3,190 1,621

DOI: 10.1371/journal.pbio.0030314.t001

studies, confirmed in this analysis, have shown that this is not
generally true for vertebrate multigene families [12]. As can
be seen in Figure S1, there is no peak at four for gene family
membership for any vertebrate. In fact, even gene duplica-
tions do not predominate; for each vertebrate species
considered individually, one member per cluster is the

largest category, accounting for 55%, 57%, and 59% of the
fugu, mouse, and human genes, respectively, with 53.4% of
the gene clusters having no duplication events whatsoever.
Thus, there is no signal of 2R remaining in gene family
membership, despite early anecdotal observations to the
contrary.

Figure 2. Overview of the Building of a Gene Cluster and Phylogenetic
Tree Shown by a Hypothetical Example
(A) Each circle represents a gene, labeled with the source genome
according to the first letter of each taxon—C, M, H, and F for Ciona,
mouse, human, and fugu, respectively—and further differentiated by
numeral. BLASTP was first used to search all vertebrate genes for the one
most similar to Ciona’s C1 gene, in this case the mouse gene M1. Then
other genes are recruited to the cluster if they have a higher similarity
score to M1 that that between C1 and M1, indicated here by the red lines.
The six genes shown on the right side of the diagram have some
sequence similarity to those in the cluster, but less than that between C1
and M1, so are not included. Because the vertebrates are more closely
related to each other than any is to Ciona, each cluster will include those
genes descended from a single gene in the common chordate ancestor,
having arisen by either lineage splitting or gene duplication specific to
one or more vertebrates. (See Materials and Methods for more details.)
(B) Evolutionary tree of the genes in this cluster show separate
duplications for fugu and for human. Because the maximum likelihood
method does not rely solely on sequence similarity, there is no
significance to the mouse gene being most similar to C1. The mouse
genome simply contained the most slowly evolving vertebrate gene in
this multigene family; this can be from any vertebrate taxon with
approximately equal likelihood.
DOI: 10.1371/journal.pbio.0030314.g002

Figure 3. Hypothetical Phylogenetic Tree Showing All Possible Types of
Gene Relationships and How They Are Most Parsimoniously Interpreted
Interior nodes are designated in lower case for those that simply result
from lineage splitting and in upper case for gene duplications within a
lineage. Although not shown, nodes are still scored if there is gene loss.
Phylogenetic trees for each gene family can be viewed at http://phigs.
org/, also providing a valuable tool for improving the inference of gene
function. DBFTS, duplication before fish–tetrapod split; DBPRS, duplica-
tion before primate–rodent split; FD, fugu duplication; fts, fish–tetrapod
split; HD, human duplication; MD, mouse duplication; prs, primate–
rodent split.
DOI: 10.1371/journal.pbio.0030314.g003

PLoS Biology | www.plosbiology.org 1702 October 2005 | Volume 3 | Issue 10 | e314

Vertebrate Genome Duplication

Figure 4. Phylogenetic Analysis of the Four Chordates with Drosophila as
an Outgroup
This phylogenetic tree is based on 766 concatenated single copy protein
sequences totaling 313,797 amino acid positions with branches propor-
tional to the amount of change. Numerals in bold above the branches
indicate the number of gene duplications occurring in each lineage;
numerals below indicate branch lengths.
DOI: 10.1371/journal.pbio.0030314.g004

Determination of Concordantly Duplicated Regions
To test the extent to which the 3,753 early duplication

events that are timed to the base of Vertebrata were
generated as part of larger scale, multigene duplications, we
examined the relative positions of these resulting paralogs in
the human genome (which is currently the best assembled and
annotated vertebrate genome). These results are shown in
Figure 5 (and more comprehensively in Figure S2) in which
the linear array of genes for each chromosome is used to
query for paralogs generated by any duplication event prior
to the fish–tetrapod split. It is apparent from these figures
that there is a large-scale pattern of genome segments that
are concordant in having similar arrangements of paralogous
genes. We quantified this by identifying all cases in which two
or more different early-duplicating genes are within a 100-
gene window, then, for each, querying all other places in the
genome, using a sliding window to count the number of cases
in which their respective paralogs are within both 50 genes
upstream as well as 50 genes downstream from that point.
There is a distinct pattern of having multiple chromosomes
matching with long linear stretches of paralogous genes. This
indicates that these duplications occurred in very large
segments, consistent with the hypothesis of WGD(s). Having
matches to three other chromosomal segments is the
dominant category, as can be seen by the darker coloring in
Figures 5 and S2 and in the histogram of Figure 6. These
patterns, with each genomic region corresponding in gene
arrangement to sets of paralogs in three other genomic
segments, are strong support for the 2R hypothesis.

Although the 4-fold (i.e., including the query segment)
category is the most prevalent, it accounts for only 25% of
the genome. Nonetheless, it is striking that this remains the
largest category despite approximately 450 MY of evolution.
This constitutes a strong signal of 2R, and could not
reasonably have been generated by a series of smaller
duplication events. For the latter to have generated this
pattern, multiple duplications of the same region (or its
resulting duplicates) would have to have occurred three
times, and have done so for many regions throughout the
genome. We would expect, rather, that independent, random
duplications would follow a Poisson distribution; this
contrasting situation is seen when the same analysis is done
with all human gene paralogs generated by duplication after

the split of fish and tetrapods (not shown). Even if we were to
consider the alternative of a single WGD followed by
subsequent independent duplications of large segments, it
would be difficult to explain why these would have been
predominantly 2-fold for previously duplicated regions. The
most parsimonious explanation for the observed pattern can
only be 2R.

Tetra-Paralogon Detection
To further establish 2R, we evaluated these sets of paralogs

for whether this 4-fold matching indicates that they fall into
tetra-paralogons, as illustrated in Figure 1. We formalized
this by first identifying paralogons (paralogous genomic
segments) containing the same set of at least two duplicated
gene pairs, while allowing for a maximum of 100 undupli-
cated genes in between (similar to the approach in [10]). (The
allowance of 100 genes is arbitrary, but the results are not
critically dependent on this number, which is only used to
find the blocks of paralogous genes.) We infer that duplicated
genes in paralogons are likely to have arisen from a single
duplication involving all contained, duplicated genes, and
that the unique, intervening genes have resulted from
differential gene deletions and subsequent genome rear-
rangements.
We identified 2,953 paralogous human gene pairs that are

inferred to have resulted from 1,912 genes that duplicated
prior to the divergence of the fish and tetrapod lineages (with
some gene losses also). Of these paralogous genes, 32.4% are
still in 386 detectable paralogons comprising 772 individual
genomic segments, containing from two to 42 gene pairs
(Table S1). Of these 772 genomic segments, 454 comprise
tetra-paralogons (Figures 7A and S3, Table 2) as shown
hypothetically in Figure 1, in which overlapping sets of
paralogs fall into 4-fold groups. (Unfortunately, it was not
possible for us to evaluate the hypothesis of an additional
genome duplication unique to ray-finned fish [41,42] because
of the generally poor contiguity of the fugu draft assembly.)
In contrast, when looking at the gene pairs that arose from

a duplication event after the divergence of the fish and
mammal lineages (see Figure 4), we find only 11% are
detected in paralogons in the human genome, indicating that
these duplications have less commonly included large seg-
ments of the genome (Figure 7B). This is especially interesting
in that their relative recency would make it more likely that
any large duplications would remain detectable, reinforcing
the contrast with the large-scale structure of those earlier
duplications. By looking specifically for tandemly duplicated
genes by defining them as paralogs on the same chromosome
that are separated by fewer than 10 intervening genes, we can
recognize that 50% of these human gene pairs arose from
tandem duplication, compared with 6% for the human gene
pairs that arose before the divergence of the fish and tetrapod
lineages.

Discussion

No detectible signal of WGD exists in the analysis of gene
family membership. There is no peak at four genes per family
for any of the vertebrates (Figure S1) as might result from 2R.
Presumably this results from a great number of subsequent
gene losses that have erased this signal. Likewise, the
phylogenetic timing of the duplication events is also incon-

PLoS Biology | www.plosbiology.org 1703 October 2005 | Volume 3 | Issue 10 | e314

Vertebrate Genome Duplication

PLoS Biology | www.plosbiology.org 1704 October 2005 | Volume 3 | Issue 10 | e314

F
ig

u
re

5
.

P
lo

t
o

f
th

e
G

e
n

o
m

ic
P

o
si

ti
o

n
s

o
f

P
a

ra
lo

g
o

u
s

P
a

ir
s

o
f

H
u

m
a

n
G

e
n

e
s

th
a

t
A

ro
se

f
ro

m
D

u
p

li
c

a
ti

o
n

s
P

re
d

a
ti

n
g

t
h

e
F

is
h

–
T

e
tr

a
p

o
d

S
p

li
t

T
h

e
q

u
e

ri
e

s
sh

o
w

n
h

e
re

u
se

C
h

ro
m

o
so

m
e

s
2

, 4
, 5

, a
n

d
1

0
, a

s
in

d
ic

a
te

d
fo

r t
h

e
fo

u
r p

a
n

e
ls

. (
T

h
e

c
o

m
p

le
te

s
e

t
c

a
n

b
e

s
e

e
n

in
F

ig
u

re
S

2
.)
O

n
t

h
e

x
-a

x
is

is
e

a
c

h
c

h
ro

m
o

so
m

e
a

rr
a

n
g

e
d

fr
o

m
p

t
o

q
t

e
lo

m
e

re
s.

O
n

t
h

e

y-
a

x
is

is
e

a
c

h
o

f
th

e
2

2
h

u
m

a
n

a
u

to
so

m
e

s
p

lu
s

th
e
X

a
n

d
Y

c
h

ro
m

o
so

m
e

s.
F

o
r

e
a

c
h

q
u

e
ry

g
e

n
e
o

n
t

h
e
x

-a
x

is
, a

‘‘
h

it
’’
is

s
c

o
re

d
if

t
h

e
s

u
b

je
c

t
c

h
ro

m
o

so
m

e
c

o
n

ta
in

s
a
p

a
ra

lo
g

g
e

n
e

ra
te

d
b

y
a

g
e

n
e

d
u

p
li

c
a

ti
o

n

p
ri

o
r t

o
t

h
e

fi
sh

–
te

tr
a

p
o

d
s

p
li

t.
T

h
e

lo
w

e
r

p
o

rt
io

n
o

f e
a

c
h

p
a

n
e

l p
lo

ts
t

h
e

n
-f

o
ld

re
d

u
n

d
a

n
c

y
a

lo
n

g
t

h
e

q
u

e
ry

c
h

ro
m

o
so

m
e

a
s

d
e

fi
n

e
d

b
y
p

a
ir

s
o

f p
a

ra
lo

g
s

d
e

te
c

te
d

in
a

s
li

d
in

g
w

in
d

o
w

a
n

a
ly

si
s.

S
e

e
t

h
e

M
a

te
ri

a
l

a
n

d
M

e
th

o
d

s
se

c
ti

o
n

f
o

r
d

e
ta

il
s,

b
u

t
b

ri
e

fl
y

, f
o

r
e

v
e

ry
h

u
m

a
n

q
u

e
ry

g
e

n
e

, a
w

in
d

o
w

w
a

s
c

o
n

si
d

e
re

d
o

f
5

0
g

e
n

e
s

to
t

h
e

le
ft

a
n

d
5

0
g

e
n

e
s

to
t

h
e

r
ig

h
t,

w
it

h
a

‘‘
h

it
’’
o

b
ta

in
e

d
f

o
r

th
e

s
u

b
je

c
t

c
h

ro
m

o
so

m
e

if
it

in

c
lu

d
e

s
th

e
e

a
rl

y
-d

u
p

li
c

a
te

d
p

a
ra

lo
g

s
o

f
g

e
n

e
s

o
n

e
a

c
h

s
id

e
o

f
th

e
q

u
e

ry
. F

o
u

r-
fo

ld
(

i.e
.,

in
c

lu
d

in
g

t
h

e
q

u
e

ry
)

m
a

tc
h

in
g

, a
s

e
x

p
e

c
te

d
b

y
t

h
e

2
R

h
y

p
o

th
e

si
s,

is
h

ig
h

li
g

h
te

d
in

a
d

a
rk

e
r

sh
a

d
e

o
f

b
lu

e
.

D
O

I:
1

0
.1

3
7

1
/j

o
u

rn
a

l.p
b

io
.0

0
3

0
3

1
4

.g
0

0
5