DOI: **.****/s*****-**6-0017-y
Evaluating Neanderthal Genetics and Phylogeny
Martin B. Hebsgaard,1 Carsten Wiuf,2,3 M. Thomas P. Gilbert,1 Henrik Glenner,1 Eske Willerslev1
1
Centre for Ancient Genetics, Niels Bohr Institute and Biological Institute, University of Copenhagen, Juliane Maries vej 30,
Copenhagen DK-2100, Denmark
2
Bioinformatics Research Center, University of Aarhus, Hoegh Guldbergs Gade 10, Building 1090, Aarhus DK-8000, Denmark
3
Molecular Diagnostic Laboratory, Aarhus University Hospital, Brendstrupgaardsvej 100, Aarhus DK-8200, Denmark
Received: 25 January 2006 / Accepted: 29 August 2006 [Reviewing Editor: Dr. Martin Kreitman]
Abstract. The retrieval of Neanderthal (Homo Introduction
neanderthalsensis) mitochondrial DNA is thought to
Ancient DNA (aDNA) studies have su ered much
be among the most signi cant ancient DNA contri-
criticism since they began about 20 years ago. The
butions to date, allowing con icting hypotheses on
modern human (Homo sapiens) evolution to be tested eld is still recovering from the e ects of early spec-
tacular and erroneous claims, such as that of DNA
directly. Recently, however, both the authenticity
being preserved in plant fossils, dinosaur bones, and
of the Neanderthal sequences and their phylogenetic
amber for many millions of years (for recent reviews
position outside contemporary human diversity have
see Hebsgaard et al. 2005; Willerslev and Cooper
been questioned. Using Bayesian inference and the
2005). Unfortunately, unreplicated results of sur-
largest dataset to date, we nd strong support for a
prising age continue to be published, including those
monophyletic Neanderthal clade outside the diversity
from old human remains (e.g., Adcock et al. 2001),
of contemporary humans, in agreement with the
expectations of the Out-of-Africa replacement model microorganisms (e.g., Cano and Borucki 1995;
Vreeland et al. 2000; Fish et al. 2002), and plant
of modern human origin. From average pairwise se-
fossils (Kim et al. 2004). These studies have routinely
quence di erences, we obtain support for claims that
underestimated the extent to which aDNA research is
the rst published Neanderthal sequence may include
confounded by contamination with modern DNA,
errors due to postmortem damage in the template
and are widely thought to result from such contam-
molecules for PCR. In contrast, we nd that recent
ination (Willerslev et al. 2004a; Hebsgaard et al.
results implying that the Neanderthal sequences are
2005). In recent years, a greater understanding of
products of PCR artifacts are not well supported,
su ering from inadequate experimental design and a postmortem damage and contamination has provided
presumably high percentage (>68%) of chimeric a more robust foundation for the eld, although the
authentication of studies of human remains and
sequences due to jumping PCR events.
microbes remains highly problematic (e.g., Willerslev
et al. 2004b; Gilbert et al. 2005a; Hebsgaard et al.
Key words: Ancient DNA Human evolution
2005; Willerslev and Cooper 2005).
Neanderthal DNA Bayesian inference PCR
The rst report of putative Neanderthal (Homo
artifacts DNA damage
neanderthalsensis) mitochondrial DNA (mtDNA)
from the type specimen (Feldhofer I [Krings et al.
1997]) was a rare example of a remarkable aDNA re-
sult obtained using very strict criteria for authenticity,
including the independent replication of results and
tests of biochemical preservation (Cooper and Poinar
Correspondence to: Eske Willerslev; email: ***********@**.**.**
51
2001; Hofreiter et al. 2001a; Paabo et al. 2004; Wil- (Caldararo and Gabow 2000; Schmitz et al. 2002) and
lerslev and Cooper 2005). The result is convincing as might possibly result from postmortem damage
the Neanderthal sequence di ers from any known (Hansen et al. 2001). Recent results from Pusch and
modern human (Homo sapiens) and chimpanzee (Pan Bachmann (2004) suggest that the original Neander-
troglodytes) sequences but is clearly human-like. Fur- thal sequences might not represent authentic
thermore, subsequent independent retrieval of similar, Neanderthal DNA but sequence artifacts.
but not identical, mtDNA from other Neanderthal Altogether these uncertainties and claims can have
specimens strongly supports the sequence s authen- severe implications for the understanding of modern
ticity (Ovchinnikov et al. 2000; Krings et al. 2000; human evolution. This paper aims to address the
Schmitz et al. 2002; Serre et al. 2004a; Lalueza-Fox et claims and reevaluate Neanderthal genetics and
al. 2005). phylogeny in an up-to-date framework. We evaluate
The retrieval of Neanderthal sequences enables the the rst published Neanderthal sequence (Feldhofer I
possibility of addressing the long-running debate [Krings et al. 1997]) with respect to damage-based
about modern human origin, something that had errors, as discussed by Gutierrez et al. (2002), and
remained unsolved in paleontological and modern further investigate the study by Pusch and Bachmann
genetic studies (Wolpo 1989; Templeton 1992). (2004) where it is claimed that the Neanderthal se-
Neanderthals have been suggested either to be (i) quences might not represent authentic Neanderthal
direct ancestors of modern man or to have contrib- DNA but sequence artifacts. The phylogenetic posi-
uted to the gene pool of today s humans (multire- tion of the Neanderthals mtDNA sequences relative
gional model [e.g. Wolpo et al. 1984; Templeton to contemporary human mtDNA is analyzed using
2002]) or (ii) to have been replaced by anatomically Bayesian inference and a newly compiled dataset
modern humans without leaving any genetic trace in together with the datasets used by Gutierrez et al.
contemporary populations (Out-of-Africa replace- (2002).
ment model [Stringer and Andrews 1988; Harvati et
al. 2003]). For a recent review of the current evidence
Materials and Methods
pertaining to this debate, see Finlayson (2005). Most
published phylogenetic analyses suggest that Nean-
Assessing Damage and Sequence Artifacts
derthal mtDNA is positioned outside the genetic
diversity of contemporary humans. This points to the
To investigate the problem of damage-based errors in the HVR1
Out-of-Africa replacement model (Krings et al. 1997, Feldhofer I consensus sequence we randomly simulated 100,000
1999, 2000; Ovchinnikov et al. 2000; Schmitz et al. sequences with the same base composition as observed in the
2002; Knight 2003), though one cannot, with the 11 positions that are variable in Neanderthal HVR1 sequences
(Table 1). This approach assumes that the Neanderthal population
limited Neanderthal sequences that are available,
is fairly homogeneous and not too genetically structured. Each of
exclude other scenarios yet (e.g., Nordborg 1998).
the 11 bases was drawn independently of each other and the
However, there are two main problems associated empirical distribution, D, of the average pairwise di erence (APD)
with the studies: (i) the use of limited contemporary to all humans was computed. The APDs between the four Nean-
human mtDNA sequences (as few as 10 [Schmitz derthal sequences and all human sequences were calculated
(Table 2), and a test was applied to determine whether the average
et al. 2002]). Such restricted sequence sampling has
obtained for the Feldhofer I sequence was extreme in D. The true
been shown to a ect the phylogenetic position of
variance of the APD is likely to be underestimated in D (but not the
ancient human mtDNA sequences (Cooper et al. mean value), because sites are drawn independently of each other.
2001); (ii) the use of analytical methods (i.e., neigh- A variance correction was therefore performed assuming the APD
bor-joining or maximum parsimony) that are, to is binomial with mean 11q as in D (the observed mean of D was
4.56, yielding q=0.414) and variance 11q(1 q) = 2.67 with
some extent, unable to account for the extreme
q=0.414 (Fig. 1). This approach is justi ed because the human
among-site variation in substitution rates and the
sequences (from large HVR1) consist of one main type comprising
levels of parallel mutations (homoplasies) that exist in 1584 of 1905 sequences, and the other types are very similar to
the human control region (Krings et al. 1997, 1999, the most frequent one. It was tested whether the APD between
2000; Ovchinnikov et al. 2000; Gutierrez et al. 2002). Feldhofer I and humans was extreme in a binomial distribution
with q=0.414 (Fig. 1).
Other important issues a ecting investigations of
To test for chimeric sequences caused by jumping PCR events
Neanderthal genetics are recent claims that the
(Paabo et al. 1989) in the datasets of Push and Bachmann (2004)
Neanderthal sequences are erroneous or simply se- and Krings et al. (1997), we used the approach of Gilbert et al.
quence artifacts. Based on disagreements in genetic (2003b), examining the clone sequences for incompatible miscoding
distances of the Neanderthal sequences to contem- lesion-derived base substitutions.
The phylogenetic analyses of the interim consensus sequences
porary humans and the age of the fossils, Gutierrez et
(ICS) from Pusch and Bachmann (2004) and the Neanderthal se-
al. (2002) argue that the rst published Neanderthal
quences were analyzed using MrBayes (Huelsenbeck and Ronquist
sequence (Feldhofer I [Krings et al. 1997]) is errone- 2003). The Markov chain Monte Carlo analysis in MrBayes was
ous. This interpretation is supported because some of run for 2 million generations with four chains, three times inde-
the positions in the Feldhofer I sequence are unique pendently. Trees were sampled every 100 generations and a 50%
52
Table 1. Variable positions among the Neanderthal mtDNA HVR1 sequences
MtDNA position
1 1 1 1 1 1 1 1 1 1 1
6 6 6 6 6 6 6 6 6 6 6
0 0 0 1 1 1 1 1 1 1 2
7 8 9 0 0 1 1 5 5 8 5
8 6 3 7 8 1 2 4 6 2 8
Sequence Frequency
Human
1 A T T C C C C T G A A 1584
2 . . C . . . . . . . . 150
3 . . . . . T . . . . . 76
4 . . . . . . . . . C . 38
5 . C . . . . . . . . . 30
6 . . . . . A . . . . . 9
7 . . . . . . . . . . G 4
8 . . . . . . . C . . . 4
9 . . . T . . . . . . . 2
10 . . . . . . . . . . C 2
11 . . . . T . . . . . . 2
12 G . . . . . . . . . . 1
13 . . C . . . . . A . . 1
14 . . G . . . . . . . . 1
15 . . C . . . . . . C . 1
Neanderthal
AF011222 G T C T T T T C G A G
AF254446 A C T C C C C T A C A
AF282971 . . T C C C C . . C .
AY149291 . . T C C C C . . C A
Note. The 11 variable positions among the Neanderthal HVR1 sequences are shown together with the corresponding substitutions in
contemporary human sequences and their frequencies.
Table 2. Mean distances between sequences of Neanderthals and contemporary humans
Average pairwise
Accession No. Name/country Ref. Age (B.P.) di erence
AF011222 Feldhofer I/Germany Krings et al. 1997 40,000 7.9092
AF254446 Mezmaiskaya/Russia Ovchinnikov et al. 2000 29,195 3.0961
AF282971 Vindija/Croatia Krings et al. 2000 42,000 4.1171
AY149291 Feldhofer II/Germany Schmitz et al. 2002 40,000 3.1234
Note. Names, countries, and ages of the Neanderthal recoveries are shown. The Neanderthal HVR1 sequences average pairwise di erences
from contemporary human sequences are also indicated. The calculations are based on the 11 variable positions listed in Table 1.
majority-rule consensus tree was produced from the last 1000 trees. Neanderthal HVR1 sequence (AY149291) not used by Gutierrez
Stationarity was checked using the command sump. et al. (2002) and two HVR1 Cro-Magnon sequences (AY283027,
AY283028) from Caramelli et al. (2003). The recently published
HVR1 sequence from Vindija (Vi-80 [Serre et al. 2004b]) was not
Phylogenetic Inference included in the analyses, as it is identical to the rst Vindija
Neanderthal sequence (Krings et al. 2000) and could be derived
To reevaluate the Neanderthal phylogeny, we took two ap- from the same individual.
proaches. In one set of analyses we used the aligned mtDNA In another set of analyses we created a dataset with 519 HVR1
control region sequence datasets from Gutierrez et al. (2002). The and HVR2 sequences of 859 aligned positions from the HvrBase
data contained a large HVR1 dataset of the hypervariable region (Handt et al. 1998; http://www.hvrbase.org). The dataset consisted
1 from the mtDNA control region, consisting of 422 aligned of 7 Chimpanzee sequences used as outgroup, the Vindija Nean-
positions consisting from 1905 contemporary human and derthal HVR1 and HVR2 sequences (AF282971, AF282972
3 Neanderthal sequences (AF011222, AF254446, AF282971). [Krings et al. 2000]), and 511 contemporary human sequences (see
Furthermore, they also included a smaller HVR12 dataset (com- Supplementary Material for the complete dataset). The Feldhofer 1
bined HVR1 and HVR2 mtDNA control region) of 843 aligned HVR sequences (AF011222, AF142095 [Krings et al. 1997, 1999])
positions, consisting of 377 contemporary human and 2 Nean- were not used to avoid bias from possible sequence errors. The 511
derthal sequences (AF011222, AF142095 and AF282971, human sequences were composed of 52 Africans, 162 Asians, 21
AF282972). Additionally, in our reanalysis, we used an additional Oceanic/Australians, 145 Europeans, and 131 Americans. The
53
invariable sites with the program PAUP* version 4 beta 10
(Swo ord 1998).
Results
Errors in the Feldhofer I Sequence
We nd that the probability of a randomly generated
sequence having a higher average pairwise di erence
than the Feldhofer I is only 1.4%. Applying the var-
iance correction we found that the probability of
obtaining an average pairwise di erence of 8 or more
(the observed pairwise di erence between Feldhofer I
and humans is 7.91; Table 2) is 3.7%. Thus, the Fel-
Fig. 1. The empirical distribution of the average pairwise di er-
dhofer I HVR1 sequence is extreme in base compo-
ence between a randomly generated sequence (see Materials and
sition compared to the Neanderthal HVR1 sequences
Methods) and the sample of human HVR1 sequences in the large
HVR1 dataset. The average pairwise di erences are mainly deter- and it is likely that the sequence is erroneous.
mined by the distance to the most common human sequence (1584
of 1905 sequences).
Arti cial Neanderthal DNA
dataset was created from the approximately 4000 taxa containing
A BLAST search (Altschul et al. 1997) demonstrates
both HVR1 and HVR2 in the database by eliminating identical or
that seven of the di erent ICS (I V, IX, and XI)
nearly identical sequences until the dataset was reduced to 511
human sequences. We used this procedure because the Markov obtained in the experiment of Pusch and Bachmann
chain Monte Carlo method that is used in MrBayes (see below) to
(2004) show 100% matches with human mtDNA
approximate the posterior probability distribution relies on con-
GenBank sequences (accession numbers AF285377,
vergence of the log likelihood and other model parameters. Con-
AF285367, AY426291, AB059953, AF519867,
vergence can be di cult to achieve when dealing with a large
AY314618, and AY314618), implying that a variety of
number of sequences and a relatively low number of unique site
patterns (Huelsenbeck et al. 2002). Additionally, we aimed to contemporary human contaminant sequences was
include the maximum amount of sequence divergence in the dataset
ampli ed. A higher frequency of transitions than
in order to be able to run the dataset within a reasonable
transversions (164/72 = 2.3), combined with the
timeframe.
higher frequency of type 2 (cytosine thymine and
To account for the large amount of parallel evolution and rate
guanine adenine, i.e., CG TA; total =95) than
variation within the HVR regions (Tamura and Nei 1993), we used
the general time reversible model of nucleotide substitution (GTR type 1 (adenine guanine and thymine cyto-
[Tavare 1986; Rodriguez et al. 1990]) with gamma-distributed rates
sine, i.e., AT GC; total = 68) mutations observed
among sites with a correction for invariable sites (GTR+G+I).
among the clone products of Pusch and Bachman
The number of gamma categories was as standard set to 4. To
(2004), is consistent with observations on postmortem
investigate the phylogenetic signals within the two HVR12 dataset
damage-derived miscoding lesions (Hansen et al.
from Gutierrez et al. (2002), we partitioned it into the HVR1 and
HVR2 sections. 2001; Hofreiter et al. 2001b; Gilbert et al. 2003a,b;
The phylogenetic analyses of the datasets from Gutierrez et al.
Binladen et al. 2006), suggesting that postmortem
(2002) were performed with the parallel version of MrBayes
damage might be involved in their data. Furthermore,
version 3 beta 4 (Huelsenbeck and Ronquist 2003) on the Bio-
at least 24 of the 35 ICS (68.6%) can be identi ed as
Cluster at the Zoological Museum, University of Copenhagen
chimeras caused by PCR jumping events according to
(http://www.zmuc.dk/). The HVR12 dataset and individual HVR1
and HVR2 partitions from Gutierrez et al. (2002) were each run the method of Gilbert et al. (2003b) in some cases,
for 15 million generations, and the large HVR1 dataset from
up to four jumping PCR events per sequence (Fig. 2).
Gutierrez et al. (2002) was run for 30 million generations. Sta-
In comparison, only 4 of 167 clone sequences (2.4%)
tionarity for these analyses was checked using the sump com-
used to generate the rst published Neanderthal
mand in MrBayes. The HVR12 dataset created for this study was
HVR1 sequence were found to contain similar evi-
analyzed using MrBayes version 3.1.1 (Huelsenbeck and Ronquist
2003) and was run for 50 million generations. Trees were sampled dence of jumping PCR events (clones A2.10, B11.4,
every 1000 generations, with a 50% majority-rule consensus tree
B11.8, and B14.9 [Krings et al. 1997]).
computed from samples after stationarity had been reached.
Figure 3 demonstrates that the diverse set of
Stationarity and e ective sample size (the number of e ectively
clones (35 ICS) obtained in the experiment of Pusch
independent draws from the posterior distribution that is sampled
and Bachmann (2004), including the XXVI se-
from) were checked with the program Tracer 1.2.1 (Rambaut and
Drummond 2004). quences, is phylogenetically more similar to the CRS
To investigate the phylogenetic signal under the neighbor-
and to each other than to any of the published
joining method (as in Gutierrez et al. 2002), we analyzed the
Neanderthal sequences. The separate position of the
individual HVR1 and HVR2 parts from the HVR12 dataset with
Neanderthals and the ICS is highly supported with a
neighbor joining using the TN93 (Tamura and Nei 1993) model
posterior probability of 100%.
with gamma-distributed rates among sites and correction for
54
mtDNA pos 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
7 8 8 8 8 9 9 9 9 9 1 1 1 1 2 2 2 2 2 3 3 3 4 4 4 5 5 5 5 6 6 6 6 6
2 2 3 8 9 1 2 3 - 6 1 3 6 7 1 2 - 3 4 0 4 6 4 5 7 3 5 6 8 0 1 2 - 9
Clone sequences
CRS T A A C T C C C - G C GA T C C - C T A C C G C A A G C A C C C - A
I . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . .
T
II . . . . . . . . . . . .. . . . . . C . . . . . . . . . . . . . . .
III . . . . . . . . . . . .. . . . . T . . . . . . . . . . . . . . . .
IV . . . . . . . . . . . .G . . . . . . . . . . . . . . . . . . . .
T
V . . . . . . . . . . . . . . . . . .. . . . ... ... . . . .
T A
VI* . . . . . . . . . . . . . . . . . .. . . ... . ... . . . .
T A
VII* . . . . . . . . . . . . . . . .. . . . ... . . .T . . .
C T T T
VIII . . C . . - - . . . . . . . . . . .. . . . ... . ... . . . .
C
IX* . . C . . . . . . . . . . . . . .. . . . ... ... . . . .
C T A
X* . . C . . . - . . . . . . . . . .. . . . ... ... . . . .
C T A
XI* . . C . . . . C . . . . . . . . .. . . . ... ... . . . .
C T A
XII* . . C . . . . . . . . . . . . C .. . . . ... ... . . . .
C T A
XIII* . C C . . . . . . . . . . . . . .. . . . ... ... . . . .
C T A
XIV . . C . . . . . . . T . . C . C. . . . ... . ... . . . .
C A C
XV . . C . . . . . . . T . . C . .. . . . ... . ... . . . .
C A C
XVI* . . C . . . . . . . T . . C . .. . . . ... . .G. . . . .
C A C
XVII* . . C . . . . . . T . . C . .. . . . ... . .G. . . . .
T C A C
XVIII* . . C . . . - . . . T . . C . .. . . . ... . .G. . . . .
C A C
XIX* . . C . . . . . . . T . . C .G . ... . .G. . . . .
C A C T T A
XX* . . C . . . . . . T . . C .G . ... . .G. . . . .
T C A C T T A
XXI* . . C . . . . . . T . . C .G . ... . .G. . . . .
C C A C T T A
XXII* . . C . . . - . . . T . . C .G . ... . .G. . . . .
C A C T T A
XXIII* . . C . . . . . T . . C .G . ... . .G. . . . .
T C C A C T T A
Fig. 2. Detection of Jumping PCR in the clone
XXIV* . . C . . . . . . . T . . C .G . .G. . .G. . . . .
C A C T T A
sequences from Pusch and Bachmann (2004).
XXV* . . C . . . . . . T . . C .G . TG . . .G. . . . .
T C A C T T A
XXVI* . . C . . . . . T . . C .G . TG . . .G. . . . .
T C C A C T T A Transitions derived from mitochondrial light strand
XXVII* . . C . . . . . T . . C .G . T. . . .G. - .
T C C A C T T A T A modi cations. Transitions derived from
XXVIII* . . C . . . . . T . . C .G . TG . . AG . . .
T C C A C T T A T A mitochondrial Light strand modi cations have been
sls_a . . C . . . . . . . T . . . . .. . . . NNN N NNN N N N N
C A C
underlined, those from Heavy strand modi cations,
sls_b* . . C . . . . . . . T . . . . .. . . NNN N NNN N N N N
C A C T
highlighted after (Gilbert et al. 2003b). Sequences
sls_c . . C . . . . . . . T . . . . .. . . C NNN N NNN N N N N
C A C
containing both Light- and Heavy-strand derived
sls_d . . C . - - - . . . T . . . . .. . . . NNN N NNN N N N N
C A C
transitions must have originated through
sls_e . . C . . - - . . T . . . . .. . . . NNN N NNN N N N N
C C A C
recombination during the PCR reaction (marked
srs_a* N N . . . . . . . T . . C .G . ... . .G. . . . .
T C A C T T A
Srs_b* N N . . . . . . T . .G . ... . .G. . . . .
T C T A C T C T T A with an asterisk).
combined HVR12, the small HVR1, and the small
Neanderthal Phylogeny
HVR2 datasets (Figs. 6a and b).
The Bayesian analyses of the four datasets from
Gutierrez et al. (2002) show that Neanderthal se-
quences are separated from modern human sequences Discussion
with a posterior probability of 100% (Fig. 4). A
Recent phylogenetic and population genetic research
schematic representation of the resulting trees for the
suggests that any genetic interchange between
large HVR1 and the smaller HVR12, and HVR1 and
Neanderthals and anatomically modern humans was
HVR2 partitions is given in Figs. 4a d. The large
very limited during the approximately 10,000 years
HVR1 (Fig. 4a), the HVR12 (Fig. 4b), and the small
(10 kyr) they potentially co-occupied the same areas
HVR1 (Fig. 4c) datasets support Neanderthal
of Europe and Asia (Currat and Exco er 2004; Serre
monophyly with a posterior probability of 100%. The
et al. 2004b) and that the Neanderthals have not
small HVR2 (Fig. 4d) datasets support Neanderthal
monophyly with a posterior probability of 63%. contributed to the mtDNA genetic diversity found in
present-day humans (Krings et al. 1997, 1999, 2000;
Analyzing the dataset created for this study shows
Ovchinnikov et al. 2000; Schmitz et al. 2002; Knight
that the Vindija Neanderthal HVR1 and HVR2 se-
2003). These issues are central to the two main the-
quences are positioned as a sister group to the con-
ories of modern human origins: the Out-of-Africa
temporary humans with a posterior probability of
replacement model, where modern humans rapidly
100% (Fig. 5). Further, this dataset shows that six
replaced archaic forms (e.g., Neanderthals) as they
sequences of African origin form a sister group to the
began to spread from Africa through Eurasia and the
rest of the contemporary humans. The neighbor-
rest of the world sometime around 100,000 years ago
joining analyses of the datasets of Gutierrez et al.
(Stringer and Andrews 1988; Harvati et al. 2003); and
(2002) con rm that the Neanderthal sequences fall
the multiregional model, where genetic exchange or
outside the sequences of contemporary humans for
even continuity exists between archaic and modern
the large HVR1 dataset. However, the Neanderthal
humans (e.g., Wolpo et al. 1984; Templeton 2002).
sequences fall within modern human variation for the
55
addition, it has been di cult to replicate the entire
Neanderthal HVR1 sequences in independent labo-
ratories (Krings et al. 1997; Ovchinnikov et al. 2000),
suggesting that preservation of Neanderthal fossils is
at the edge of what is required for successful DNA
studies. It is therefore possible that some of the
published Neanderthal DNA sequences might con-
tain errors due to miscoding lesions (Hansen et al.
2001). This type of DNA damage is of particular
concern if ampli cations start from few template
molecules, which appears to be the case at least in the
rst published Neanderthal study (the Feldhofer I
HVR1 sequence [Krings et al. 1997]).
In support of errors in the Feldhofer I HVR1 se-
quence it has been argued that the most recent Nean-
derthal specimen (Mezmaiskaya, 29 kyr old) shows a
shorter genetic distance to contemporary humans than
Feldhofer I (which is believed to be the oldest of the
Neanderthal specimens) (Gutierrez et al. 2002). How-
ever, the validity of this argument is questionable, as
the Feldhofer I fossil has recently been redated to 40
kyr (Schmitz et al. 2002), and the young age of the
Mezmaiskaya fossil is debated (Skinner et al. 2005).
Additionally, it has been noted that the Feldhofer I
HVR1 sequence harbors four unique substitutions
Fig. 3. The phylogenetic position of the Neanderthal sequences (positions 107, 108, 111, and 112) possibly due to
relative to the clone consensus sequences generated by Pusch and
postmortem damage accumulated during ampli ca-
Bachmann (2004). The XXVI clone consensus sequence is the one
tion (Caldararo and Gabow 2000; Schmitz et al. 2002;
claimed to be Neanderthal-like.
Hansen et al. 2001). Using a maximum damage-based
error rate of 0.06%, Hofreiter et al. (2001b) reject
In this paper we have investigated the genetic major errors in the Feldhofer I sequence. However, the
a nities of the Neanderthals to anatomically modern rate might be underestimated because they do not take
humans. First, we have evaluated whether the rst into account the possible presence of damage hotspots
published Neanderthal sequence (Feldhofer I) is in the human D-loop (Gilbert et al. 2003a) and the
erroneous (Gutierrez et al. 2002). Second, we have error rate is calculated from the consensus of only three
investigated whether the Neanderthal sequences are Neanderthal sequences.
sequence artifacts (Pusch and Bachmann 2004). Comparisons of the average APD between the
Finally, with our re ections on the rst two questions Feldhofer I sequence and the sequences of contem-
in mind, we have readdressed the controversial ques- porary humans with the APD of randomly generated
tion about the phylogenetic position of the Neander- Neanderthal sequences indicate that the Feldhofer
thals. I HVR1 sequence is extreme in its genetic composi-
tion. It is therefore likely that the Feldhofer I HVR1
Errors in the Feldhofer I Sequence sequence is erroneous and we cannot exclude that at
least this sequence is modi ed due to postmortem
One explanation for the unresolved position of damage (see Errors in the Feldhofer I Sequence,
Neanderthals among anatomically modern humans is under Results).
that the sequence data might be considered unreliable
due to the degraded nature of the Neanderthal Arti cial Neanderthal DNA
specimens and their DNA (Gutierrez et al. 2002).
Biochemical analyses for investigating the preserva- Instead of the Neanderthal sequences being a ected
tion condition of excavated Neanderthal bones and by damage, a recent study suggests that their unique
teeth indicate that most of the specimens are unlikely substitution patterns are caused by PCR artifacts.
to yield any endogenous DNA (Serre et al. 2004b). Pusch and Bachmann (2004) report that 35 di erent
The majority of samples that have yielded putative mitochondrial HVR1 sequences, including a group
Neanderthal DNA have only enabled PCR ampli - containing 7 substitutions that are in combination
cation of mtDNA in the 50-base pair (bp) size range characteristic for the Neanderthals (i.e., clone XXVI;
(Serre et al. 2004b; Lalueza-Fox et al. 2005). In Fig. 2) can be ampli ed from a single sequence of
56
Fig. 5. Majority rule consensus trees from the MrBayes analyses
showing the phylogenetic relationship among the seven Chimpan-
zee sequences used as outgroup, the Vindija Neanderthal (Krings
et al. 2000), and 511 contemporary humans.
Fig. 4. A schematic representation of the majority rule consensus
trees from the Bayesian analyses of the datasets from Gutierrez
et al. (2002) showing the phylogenetic relationship between the
Neanderthals and contemporary humans. (a) The large HVR1
dataset, (b) the HVR12 dataset, (c) the HVR1 partition of the
HVR12 dataset, and (d) the HVR2 partition of the HVR12 dataset.
modern human mitochondrial DNA (matching the
Cambridge Reference Sequence; CRS [Anderson
et al. 1981]) if the PCR reaction is spiked prior to
ampli cation, with 14 di erent aDNA extracts of
non-Neanderthal origin. The authors thereby indi-
rectly imply that the published Neanderthal se-
quences could be explained in this manner and may,
in fact, not represent authentic Neanderthal DNA.
Fig. 6. Neighbor-joining tree using the datasets from Gutierrez
However, as shown in Fig. 3, the diverse set of clones et al. (2002). (a) The large HVR1 dataset and (b) the combined
HVR12 and the individual small HVR1 and HVR2 partitions from
(35 ICS) obtained in the experiment of Pusch and
the HVR12 dataset.
Bachmann (2004), including the Neanderthal-like
XXVI sequences, is phylogenetically more similar to
the reference sequence and to each other than to any
mtDNA GenBank sequences. This strongly suggests
of the published Neanderthal sequences. The separate
that a variety of human contaminants is ampli ed in
position of the Neanderthals and the ICS is highly
supported, with a posterior probability of 100%, and the experiment. Furthermore, the higher frequency of
transitions than transversions combined with the
all of the arti cially generated sequences are therefore
higher frequency of type 2 than type 1 mutations (see
clearly distinguishable from the published Neander-
Materials and Methods) among the clone products of
thal sequences.
Pusch and Bachman (2004) is consistent with the
Another interesting issue is that regular BLAST
presence of damage-based misincorporation in the
searches (Altschul et al. 1997) reveal that seven of the
template DNA (Hansen et al. 2001; Gilbert et al.
di erent ICSs obtained by Pusch and Bachmann
(2004) show 100% match with di erent human 2003a,b, 2005a; Willerslev et al. 2003). This could
57
Table 3. Sequence comparisons of mtDNA HVR1
MtDNa position
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
7 8 8 8 8 9 9 9 9 9 1 1 1 1 2 2 2 2 2 3 3 3 4 4 4 5 5 5 5 6 6 6 6 6
2 2 3 8 9 1 2 3 6 1 3 6 7 1 2 3 4 0 4 6 4 5 7 3 5 6 8 0 1 2 - 9
Sequences
Human
CRS T AACTCCC GCGATCC CTACCGCAAGCACCC- A
XXVI . . CTC. . . C. . ATC. . CT. GT. ATG. . . G. . . . .
Neanderthal
AF011222 A C C C C C C C A C A C T C C T T C T A A C A C A A G C C C T A
AF254446 A C C C C C C C A C A C T C C T T C T A A C A C A A A C C C T A
AF282971 A C C C C C C C A C A C T C C T T C T A A C A C A A G C C C T A
AY149291 A C C C C C C C A C A C T C C T T C T A A C A C A A A C C C T A
Human GenBank
AY244327 C C T . . . .
AF500966 . . . . A . .
AY217640 . . T . . A .
U54381 . . . . . . G
AF500983